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ABSTRACT 


The precision of a measurement of a gas-chromatographic peak 
depends on the intensity and character of the base-line noise, the 
Shape and size of the peak, the prior information, and the method of 
estimation. An estimation method is presented that approaches the 
highest precision capable of any method used with a given gas- 
chromatographic system. Unlike most other Sane ie methods, the 
method has the ability to use all the information available to 
estimate the peak parameters, including prior information as well as 


information from the chromatogram. 


The base-line noise of many detectors can be statistically 
characterized by a power-density spectrum. In such systems, the 
estimation of peak parameters is most easily carried out in the 
frequency domain on the Fourier transform of the chromatogram. 
Although any mathematical model that can completely describe a peak 
is suitable for estimation purposes, a peak model that uses cumulants, 


which are related to moments, has certain mathematical advantages. 


Digitization of a continuous chromatogram involves the loss 
of some information, but the loss can be controlled by the choice of 
the sampling interval, the sampling time, and the quantization interval. 
Essentially all the information in a chromatogram about a peak lies in 
a section of the chromatogram, centered about the mean of the peak, 
which has a length equal to eight standard deviations of the peak plus 


the autocorrelation width of the base-line noise. 
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The estimation method was verified experimentally by 
analyzing mixtures of hydrocarbons on a simple isothermal gas 
chromatograph employing a thermal conductivity detector. Standard 
deviations of 0.003 relative for peak area and 0.02 index units 
for retention index were obtained. As shown by the analysis of 
computer-generated peaks, these standard deviations could be reduced 
to 0.0005 and 0.002 if the experimental conditions were carefully 
controlled. A peak model with seven adjustable parameters adequately 
fitted moderately sized peaks of aliphatic hydrocarbons, but was 


inadequate for large peaks and peaks of aromatic compounds. 
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CHAPTER 1 


THEORETICAL 


In gas chromatography the ultimate concern is with the peak(s) 
that are representations of the distribution of solute molecules ina 
gas. Information about a peak is obtained from a chromatogram, which 
ordinarily is a voltage-time representation of the peak. The informa- 
tion in a chromatogram, however, is corrupted by random fluctuations, 
which we call base-line noise, that are not directly attributable to 
the peak. Base-line noise adversely influences the precision of 
measurement of peak parameters. This thesis is concerned with the 
application of statistical communication theory and estimation theory 
to the minimization of the effects of base-line noise on the estimation 
of peak parameters. This study is not concerned with the source of 
base-line noise, but only its effect. Moreover, this study is not 
concerned with distortions in the observed peak caused by the detector 
itself an any other determinate errors that can be corrected by 


Calibration. 


The theory of estimation, and the mathematical techniques 
used to describe peaks and noise that are presented here have been 
successfully applied in many fields outside of gas chromatography: 
economics, communications, radar, and seismology, among others. Since 
these fields are outside the normal range of contact of people working 
in gas chromatography, this chapter attempts to provide a general, 
though incomplete, discussion of these methods and a brief guide to the 


literature. Most of the general material is described more completely 
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2 
in References 1 to 4. The primary purpose of this chapter, however, is 
to describe and explain the techniques used in Chapter 2 for the 
analysis of peaks obtained from an isothermal gas-chromatographic 
system with a thermal conductivity detector. Undoubtedly these methods 


are applicable to many areas of chemistry outside of gas chromatography. 


In order to distinguish a peak from extraneous noise or non- 
peak signal components, the shape of the peak must be defined at least 
as accurately as the precision of the measurements warrant. For 
example, to define a peak as something that has area, retention time, 
and width is not enough; noise waveforms may have these attributes also. 
Since a general closed mathematical definition (or model) for the shape 
of a peak has not yet been developed, a model consisting of a series of 
terms of decreasing importance should be fruitful provided the number 
of significant parameters required is not impractically large. Series 
representations of variously shaped peaks can be developed from several 
basis functions: Gaussian, lognormal, Poisson, Beta, or Gamma (5). 

The Gram-Charlier type A series, based on the Gaussian function, has 
received the most attention in the literature on gas chromatography 
(6-9), and a special form of it (Edgeworth's form) has been chosen for 
study here. The third and higher coefficients of Edgeworth's form of 
the Gram-Charlier type A series, defined below, are called cumulants. 
Cumulants have been widely used for the characterization of statistical 
distributions. They are related to moments but they are mathematically 
easier to deal with and we believe they are more fundamental than 
moments. Although cumulants do not have direct physical significance 


in the same way that peak area is related to quantity of a solute, they 
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have theoretical significance (6-9) and may provide qualitative 


information about a solute (9,10). 


For the purposes of this study it is assumed that the 
essential characteristics of base-line noise can be described by the 
Statistical distribution of ordinate values of the points on a base 
line. Specifically, base-line noise is assumed to be a stationary 
random process with a normal or Gaussian probability-density function. 
Since the base-line noise of most gas-chromatographic detectors results 
from several independent sources, the assumption of a normal probability- 
density function should be accurate”. When the statistical behavior 
of noise does not change with time, the noise is said to be stationary. 
The assumption that base-line noise is stationary during the elution of 
a peak should hold when the experimental parameters that affect the base 
line are not changing with time, as in isothermal gas chromatography; 
but it may not hold in programmed temperature gas chromatography if the 
intensity of the base-line noise increases as the column temperature 
increases. Noise with an intensity or variance that increases with 
increasing amplitude of the signal (as during the elution of a peak) is 
nenstattonary andv1s usually difficult to treat inthe presence of 
stationary noise. Noise that depends on the amplitude of the signal may 
constitute a major part of the noise of ionization detectors and 
detectors that measure the radioactivity of a sample, but should be 
negligible for thermal conductivity detectors. 


* . . 
The central-limit theorem states that subject to certain 


very general conditions, the sum of random variables becomes more and 
more normal as the number of independent random variables increases 


(11) . 
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The fluctuations in measurements of peak parameters are not 
entirely due to base-line noise. For example, small random variations 
in column temperature and column inlet pressure can cause fluctuations 
in peak retention times and shapes without affecting the base line. In 
concentrating on the effects of base-line noise the approach taken here 
is complementary to that of Goedert and Guiochon (12). These authors 
have studied the effects of random variations in the experimental 
conditions, which are classified here as non-base-line noise, on the 
measurement of peak area and peak height. The total variance of peak- 
parameter measurements is a combination of the variance due to non- 
base-line noise and that due to base-line noise. The variance due to 
base-line noise alone can be estimated from the statistical character- 
istics of a section of base line, and this estimated variance used to 
separate the contributions of the two types of noise to the total 


observed variance of the parameters. 


Peak retention time is commonly measured from the peak maximum, 
which is only a single point from the peak profile. Points on the sides 
of the peak, however, also contain information about retention time. 

If the entire peak is considered in the calculation of peak parameters, 
a decrease in the variance of the parameter estimates should be 
attainable. Noise always adds variance to a parameter estimate, in an 
amount that depends in part on the method of estimation or measurement. 
The estimation method developed in this thesis is optimal or nearly 
optimal in the sense that the variance of the parameter estimates 
approaches the ultimate minimum variance. This is achieved by using 


all of the information in the chromatogram, not merely the information 
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obtainable from a few points such as the peak maximum and the peak 


width at the base line. 


The computation required in the estimation method is 
sufficiently complex so that an offline digital computer is necessary. 
Using a computer means that the data must be digitized. Some of the 
information contained in a continuous signal is invariably lost on 
digitization. Nevertheless, interval-area samples taken at constant 
intervals of time (13) contain almost all the information ina 
chromatogram provided the sampling interval is eee than a certain 
minimum value. Also, the information lost from quantization is 


negligible if the quantization interval is small relative to the 


intensity of base-line noise. 


The notational conventions of Bracewell (1) have been 
followed whenever possible. One of these is that the Fourier transform 
of a function is capitalized. Others are the use of the convolution 
symbol * , and the letters Ill and II from the Russian alphabet, 
which have a pictorial similarity to the functions they represent. 
Three different time and frequency scales are used. Time in seconds 
is represented by the symbol Tt , and frequency in the corresponding 
scale (hertz) by v . Time in peak-width units and in sampling- 
interval units is represented by the symbols ©¢ and 2 , and frequency 


in the corresponding scales is represented by n and s. 
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PEAK MODELS 


Peak shapes can be most easily compared when peaks are 
reduced to unit area and unit width, and are centered at the origin. 


Defining the area, A , retention time, ¢t. , and width, o , as follows: 


O 
EE GG cler (1) 
f foe} 
eri J tufts) at (2) 
ae a 2 
ei Nein (eos) flt-t) dt (3) 


a dimensionless peak-shape function, y(t) , can be derived 


from a peak on a noise-free chromatogram, f(t) , through the formula 
y(t) = f(t)o/A (4) 


where © is dimensionless reduced time: 


ee (t-t )/o (5) 


The above definitions of the retention time and width of a peak are the 
mean and the standard deviation. These are different from the more 


common definitions as the time of the maximum and the width at the base 


or a certain fraction of the height. The peak models discussed below 


are simplified when the above definition of a peak-shape function is 
used, and this is the primary reason why we define retention time and 


width as the mean and standard deviation. 
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7 
With a noise model applicable to most detectors (described 

in the next section) the parameter estimates are best calculated on the 
Fourier transform of the chromatogram. Moreover, the mathematical 
expression for the peak-shape model defined below is simpler as a 
function of frequency than as a function of time. Even though the 
Fourier transform operator may be complicated, the advances in Conan 
algorithms for numerical Fourier transforms make the use of the operator 
of little concern (14). The Fourier transform of a standardized peak 
shape function is 


foe) 


Yin) = { exp[-t2mz] y(t) az (6) 
—O0O 
where z is the square root of -l, and n is frequency ina 


dimensionless standardized scale. 


Parameters 

The three parameters area, A , retention time, bo , and 
peak width, o , along with a peak-shape function, completely and 
uniquely describe a peak. Knowing the values of the three parameters 
is enough to distinguish a particular peak from any other peak only 
when the peak-shape function is previously known. Because of the 
complexity of the gas-chromatographic process and the difficulty of 
controlling experimental conditions, accurate prior knowledge of the 
peak shape function is usually not available. Therefore, in the 
mathematical model of a peak, parameters in addeeLon. tO. Am. bo 5 mand 


Oo are required to describe the peak-shape function. For estimation 


purposes, any set of parameters that can completely and uniquely 
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8 
describe the peak-shape function will serve. It is not necessary that 
the parameters have an easily interpreted physical meaning, especially 
Prony A; by » and oO are required for subsequent calculations, but 
it is desirable. Three possible sets of parameters of a peak-shape 


function are samples from a previous chromatogram (15), moments (Se 


and cumulants. 


me a2rst, Second, thirds; ...° “moments! of a peak-shape 
function may be defined as derivatives of the Fourier transform of the 
peak-shape function: 


a) pie) 
LF (-121) ~D [Yin)] 9 (7) 


where Ds is the differential operator, Beano , and Yy is equal 
eOnwL,2,5;,00. © Ln contrast; the cumulants of a peak-shape function 
are defined as derivatives of the logarithm of the Fourier transform 
of the peak-shape function: 


K. = (-421) "D {In Y(n)] 6 (8) 
4) n= 


The moments and cumulants defined by Equations 7 and 8 are ordinary 
moments and cumulants of the Rew aneie function defined by Equation 4. 
These dimensionless numbers are referred to here as the standardized 
moments and cumulants of a peak. The first three standardized 


cumulants are equal to the first three standardized moments. But, a 


cumulant of a higher order, Yr , is a function of the lower-order 


moments from 1 to r.. Tables of the algebraic formulas relating 


moments to cumulants are given in Kendall and Stuart (16). The first 
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9 
cumulant of a standardized peak-shape function is equal to zero, and the 
second is equal to one. The third and the fourth standardized cumulants 


are sometimes called skew and excess. 


An infinite number of parameters in any one of the above sets 

will completely specify any physically possible peak-shape function. 

But an infinite set of parameters is unnecessary because in the presence 
of noise the peak-shape function can never be specified exactly. 
Therefore, if the parameters are arranged in order of decreasing 
Significance (increasing statistical variance or decreasing information 
content), only the lower-order parameters up to the number that can 
describe the peak-shape function to a required level of significance 
need to be estimated. The parameters that are estimated are called the 
adjustable parameters, and the remaining parameters in the infinite set, 
which are given an a prtori value, are called the nonadjustable 


parameters. 


The estimation of a large number of parameters is undesirable 
as well as unnecessary because, as the number of parameters to be 
estimated increases, the peak-shape function becomes increasingly 
general and thus able to fit a wider range of differently Shaped peaks. 
A model that is too general is likely to include nonpeak components such 
as noise or interfering peaks along with the true peak in a chromatogram. 
One way of introducing into the peak model the prior information 


needed to distinguish a true peak from noise is by limiting the number 


of adjustable parameters. When there is little prior information about 


the model required for a particular peak, a compromise must be made; 


the number of adjustable parameters should be enough so that the peak 
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10 
model can fit the peak, but not so many as to fit any nonpeak 


components that may be present. 


In this study cumulants are chosen to be the parameters for 
peak-shape functions. Since moments can be calculated from cumulants, 
cumulants have the same theoretical interpretation as moments, but 
they are easier to manipulate [they are additive under convolution, 
(6)]. Simplicity is the only objective reason why we choose cumulants 
over moments. The model incorporating cumulants as parameters is 
chosen with the prior information that in simple ideal chromatographic 


systems the peak-shape function approaches the Gaussian function. 


Peak Shapes 

A peak-shape function can be represented by a series of terms 
that are functions of the cumulants of the peak-shape function by use 
of the general Charlier derivative-series method (5). Given the 
Fourier transform of some basis function, 9%(n), the series expansion 


of the Fourier transform of a peak-shape function, Y(n), is 


Y(n) = o(n) expl } (Ky) (-t2mm)"/28] (9) 
y=l 
where K and Y, are the r'th cumulants of the peak-shape function 
and the basis function. The basis function should be chosen so that 
only the first few cumulants of the peaks being investigated need be 
estimated. For example, the Poisson function might be used for 
moderately tailed peaks. When a Gaussian basis-function is used, 


Equation 9 becomes 
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Y(n) = exp[ ) K (-t2mn)”/p1) (10) 
r=l1 
because all the cumulants of a Gaussian function are zero, except 
for the second, which is equal to one. An equivalent form of 
Equation 10 can be obtained by expansion of the exponential and 
collection of powers of 7n : 


foe) 


Y(n) = [1+ ) U,(-t2m)”/r!jexpl (-i2m)*/2] (11) 
p=3 


2 
where U, and U., are equal to Ky and Ky , and Ue = K. + 10K, : 


UL = Ko =e 35K kK. Fi Ue = Ky + S6K.K, + 35K4 aerate Equation 11 
defines the Gram-Charlier type A series (GCA series), and 
Equation 10 defines Edgeworth's form of the GCA series (17) Ae ere 
most frequently used series expansion of a chromatographic peak (6-9) 
is the time domain version of Equation 1l: 


° aie) 
Meee, URAC) (12) 


r=3 


where hc) is the Tchebycheff-Hermite polynomial: 


3 2 
We el ea (13). 
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Equation 9 , carried to the sixth cumulant, is used as a model in 
Chapter 2. 


Although the GCA series and the Edgeworth series are 


formally equivalent, terminating them after a finite number of terms 
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involves different prior information about the nonadjustable parameters. 
Consequently, the effects of the termination on the precision and 
accuracy of the adjustable parameters may be different. The coefficients 
of the sixth and higher terms of the GCA series are functions of 
lower-order cumulants, whereas the coefficients of all the terms of 
Edgeworth's series are the cumulants themselves. Since the high- 
order coefficients of the GCA series contain information about lower 
order cumulants, neglecting a high-order coefficient would involve the 
loss of some information about the lower-order cumulants. For example, 
suppose the third cumulant of a peak-shape function were equal to one 
and the cumulants of higher order were equal to zero. This peak- 
shape function is a possible representation of an actual peak, but since 
the peak would be partly negative, it would be unusual. The series 
representation of the peak in Edgeworth's form could be terminated after 
the term involving the third cumulant. On the other hand, the 
representation of the same peak in GCA form would have to be continued 
at least to the sixth term because from the contribution of the third 
cumulant, the sixth term would be Squal. to ten, a value that might be 
Siguiticant even 1f the sixth cumulant were insignificant: » For 
simplification of the effects of termination on parameter estimation, 
therefore, any one term of a series model of a peak-shape function should 


depend on only one parameter. 


The coefficients of the GCA series being functions of more 
than one cumulant does not mean that this series should not be used for 
the representation of peak-shape functions. For certain peak shapes in 


certain types of noise, the coefficients of the GCA series, se ; 
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could conceivably approach insignificance faster than cumulants and 


would therefore be preferred to cumulants. 
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A NOISE MODEL 


Errors in the estimates of peak parameters due to base-line 
noise are caused by the unpredictability of the noise. But, because 
all physical signals are bandlimited (18), base-line noise is never 
completely unpredictable. Given a point on a base line, a second point 
occurring some time later will likely be within a range of the first 
point that depends on the intensity and character of the noise and 
the separation of the points. By giving the most weight to points on 
a base line nearest a peak, a chromatographer uses this relation, 
sometimes unconsciously, when he examines the base line under the peak. 
When noise is stationary, so that its average properties do not change 
with time, and when the mean of the noise is zero, the randomness as 
well as the predictability, or nonrandomness, of noise is characterized 
by the autocovariance function, C(t) , which is the average of the 
product of the ordinate values, x(t) , of two points on a base line 
Taken T Units of time“apart. (19). 

- 1/2 

e (t) = lim = f (t) a(tt+t) at (14) 

i foo t ° _p/2 
where t is called the lag. Sometimes the terms autocovariance and 
autocorrelation are interchanged, but we reserve the term autocorrela- 
tion for a normalized autocovariance (see Equation 49). The auto- 
covariance at zero lag is the variance of the noise; as the lag becomes 
large, the autocovariance should approach zero, indicating that points 


far apart are independent (covariance independent). 
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Power-Density Spectrum 

The Fourier transform of an autocovariance function is the 
power-density spectrum*, Ca) +) LE, tf Janderso Navertadimensi one vot 
volts and seconds, Ce will have dimensions of volt squared per 
hertz. The power-density spectrum of noise is the variance density 
per unit frequency, hence the area of a power-density spectrum, which is 
the integral of the variance density over all frequencies, is equal to 
the variance or power of the noise, ec, (0). White noise is noise 
that has a constant or flat power-density spectrum over a wide frequency 
range. The autocovariance function of white noise consists of a 
spike-like function which rapidly approaches zero on both sides of the 


origin. 


The method for estimation of peak parameters described later 
requires the use of the probability-density function (p.d.f.) of the 
Points On a Dase Line. It the p.d.£. of the noise is normal, “as 
assumed here, then the noise is completely characterized by its power- 
density spectrum. In the limit of continuous sampling so that integrals 
Gan replace summations, the joint p.d.f. of a set of points, £ , on 
a section of base line is 


: Ixv) |? 
Pat) = b exp[-1/2 f a7 sled GF halon v] (15) 
—CO aG 


where X(v) is the Fourier transform of the section of base line being 


*n power-density spectrum may also be called a power spectrum 


or aspectral density as in Reference 3. 
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studied. The existence of the integral in Equation 15 requires that 
the power-density spectrum, C (Y) , be greater than zero for all 
frequencies. This causes no difficulty in practice when the data are 
quantized because quantization introduces white noise with a power- 
density spectrum greater than zero over the range of frequency being 
considered. The constant, b , in Equation 15 is independent of the 
particular section of base line being studied; therefore its evalua- 


tion is unnecessary. 


A point on a base line in the time domain is dependent on all 
the other points for which the autocovariance is greater than zero. 
The dependence may extend over a large number of points and over a 
wide period of time. Calculations involving stationary normal noise 
are considerably simplified in the frequency domain because the 
Fourier transform of highly interdependent base-line noise results in 
a set of effectively independent points, the variance of each point, 
X(v) , being equal to CY) . If the points on the base-line were 
effectively independent of each other, that is, if the noise were 
white over the frequency range of interest, of course, there would be 


no need to do the calculations in the frequency domain. 


The probability-density function of a section of base line, 
P, (£) , is used only formally in the estimation calculations; its 
Sve value is not determined. Nevertheless, the value of Py) 
for a particular section of base line is proportional to the probability 
that the base line is caused by noise alone. For example, if the value 


of p(x) is unusually low, the chromatogram probably contains a 


component that is not noise, possibly a peak. If a peak is 


subtracted from the base line, the value of P,, (£) might 
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become reasonably large, and the increase may indicate that the 


original chromatogram contained a peak similar to the one that was 


subtracted. 


Power-Density Spectrum Estimation 

The calculation of the integral in Equation 15 requires 
prior knowledge of the power-density spectrum. The definition of 
the autocovariance function (Equation 14) suggests that the power- 
density spectrum may be estimated from the Fourier transform of a 


section of a chromatogram that does not contain a peak (20) 
2 
Cv) Awa whl SKE (16) 


where Y is the interval of time over which data is taken. The 
section of base line used for the estimate of the power-density 
spectrum could be taken some time before or after the peak being 
examined; all that is required is that the noise be stationary enough 
so that the noise near the peak will not differ greatly from the noise 
used to estimate the power-density spectrum. The variance of the 
spectral-density estimate, CY) in Equation 16, is equal to Creye 
and is independent of the amount of data used. The precision can be 
improved by averaging the spectrum over a small range of frequency. 


The averaging process resulting in Cry) , called smoothing or 


convolution, is described by the equation: 


Crlv) = f C (ud) Wiv-u) du (17) 


—_—C 


where W(v) is a weighting function called a spectral window (20). 
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18 
The convolution operation is conveniently represented by the symbol 


* . With this symbol, Equation 16 can be written 


Cav) = Cty) x W(v) (18) 


The variance of the smoothed spectral-density estimate, Catv) eas 


co 


variance of Catv) ve Cl(y)? bee neyo ke au (19) 
After convolution with W(v) , the variance of the power-density 
spectrum decreases with increasing amounts of data and can be made as 
small as desired. But, the increase in precision obtained by smoothing 
is at the expense of accuracy; the estimate will be biased by an amount 
proportional to the curvature of Ge AY and to the width of the 
Spectral window. The bias can be decreased by filtering the data to 
give a nearly white spectrum before smoothing (21). We apply this 
ftziter to the data in the time domain to minimize end effects as 
discussed later in the section on sampling. In the time domain, 


filtering resulting in wx’(t) is represented by convolution 
iis oleae L) (20) 


where h(t) is the impulse-response function of the filter. Since the 
mean of a base line in gas chromatography is arbitrary, it has a high 
variance indicating that the power-density spectrum of base-line noise 
will be extremely high at zero frequency. To whiten this intense 


noise at zero frequency, the filter h(t) should have zero area 


(Figure 1, II). The digital filter, [-1,2,-1], gives an approximately 
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flat spectrum for the base-line noise of a thermal conductivity 
detector (see Chapter 2). The compromise between bias and precision 
that governs the choice of a spectral window is not critical when the 
spectrum is white. We have used a simple digital spectral window, 


(4575574) . 


Convolution in the time domain (Equation 20) corresponds to 
multiplication in the frequency domain; therefore, the power-density 
Spectrum can be regained from the smoothed spectrum by multiplying the 
smoothed spectrum by the inverse of the absolute square of the Fourier 
transform of the impulse-response function of the filter. At zero 
frequency both the spectrum of the filtered data and the filter- 
transfer function, #H(0) , are zero leaving the power-density spectrum, 
C, (9) , undefined. Since X(0) , the mean of the base-line noise, has 
a high variance, C60) is usually extremely large; therefore, its 


reciprocal, needed in Equation 15, is defined to be zero. 


Figure 1 shows the steps involved in the estimation of a 
power-density spectrum. First, filter the data in the time domain 
using a filter that will give a nearly white spectrum (Figure l, I to 
III). Fourier transform and calculate the power-density spectrum using 
Equation 16. Smooth according to Equation 18 (Figure 1, IV to Vi) 
and then multiply by the inverse of the square of the transfer function 
of the filter (Figure 1, VII and VIII). These operations are 
summarized by 


(|x'(v)|? « W(v)] le(v)|~? (21) 
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FIGURE l 


] I] IT] 


Base line Filter Filtered 
noise noise 


x< 


Voltage — 


Time ae 
IV V VI 


Raw Spectral Smoothed 
spectrum window spectrum 


—|>< 


re 
a” 
= 
o 
Dy 


VII VIII IX 


! I 
Filter | Final spectrum | Biased and 
correction | estimate  Aunsmoothed 


Ry \ 


7\i 


Power 


Frequency ——> 


{ mr) 
- \ 
| 
| uth 
_ e 
} 
Al 


bac Ce wes ri 


a ney 


aitnees a 


iW 8 We ts fol _ 
2 7 Tt " 
7 s { : as y es 4 a, 
; - a , > " 
} i oe 
j ; ry ' r : cu 
! 1! iF ve! 
Zl i 
: ] i : — § 7 
; he le ae wes 
; ’ iis 
Seana | 
’ 4 , 


rerya 3e a 
whan 
j 
on a "3 
| hy 
fe i’ 
a » ' 
: | 
hd 


ity 
clay ome Mei RIK SS 
avmmites 


Re eH a 


ais 7 
pe wr : 4 


i ahd) 
1 


= yanaupes. 
: < 


te 


$5 ’ 7 
, Ja a 


i 


21 


Figure 1. Estimation of power-density spectrum (schematic). 

£, section of base-line noise; II, impulse-response function of a 
whitening filter; III, filtered data resulting from the convolution 
of I and II; IV, power-density spectrum of the filtered data; 

V, a Spectral window; VI, smoothed power-density spectrum obtained 

by convolving V with IV; VII, reciprocal of the absolute square 

of the impulse-response function of the filter in II. The dashed 
line in VIII is the final power-density spectrum of the base-line 
noise obtained by multiplying VI by VII. The solid curve in VIII 
and IX is the true power-density spectrum obtained as T 

approaches infinity. The upper-dashed curve in IX is the unsmoothed 
Spectrum resulting from the multiplication of IV and VII. The 
lower-dashed curve in IX is the biased spectrum obtained by smoothing 
mitnout whitening. Since all of the curves except I and III are 
symmetric about the ordinate axis, the negative parts of the curves 


have been omitted. 
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with 


Two useful, simple characterizations of noise are the power 
or variance and the autocorrelation width. The power is defined by 


the equation 


Was if Cty) dv (22) 


and the autocorrelation width by 


T, = C0) /W, (23) 


One half the autocorrelation width is an indication of the separation 
in time between independent or uncorrelated points in the noise. For 
the calculation of ee and T. in practice, it can be assumed that 


peo) © C.(1/T). 
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DIGITIZATION 


Converting an analog signal such as the voltage from a 
detector into a form that can be used by a digital computer involves 
two distinct processes. First, a finite number of samples are taken; 
second, the samples are quantized. The resulting digitized Signal is 
a finite series of numbers each having a finite number of significant 
digits. The effect of these processes on the precision and accuracy 


of sampled data is described in detail below. 


Sampling 

The application of the theory presented in (1) is discussed 
in the following section. The information contained in a chromatogram 
and the information contained in a set of samples of the chromatogram 
can be compared when the set of samples is treated as a function of 
continuous time rather than discrete time. The delta function, 46(€), 
which is a generalized function and not a function in the usual sense, 
but can be treated as if it were, represents a single sample. The 
delta function can be defined as the limit of a Gaussian function as 
the standard deviation approaches zero; although it is infinitely high 
and infinitesimally narrow, it has unit area. Figure 2, I shows a 
delta function at time t’ as an arrow of unit height. Multiplying 
a data function, x(t) , by 6(t-t') results in a function defined 
for all time that represents a sample of the data at time t! 
(Figure 2, II). The operation of taking samples of a Function). “soe, 


Reeeonstant intervals of time, At seconds apart, is described hy 


x(a) = \ 68(a-n)a(a) (24) 
S) 
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FIGURE 2 


x(t) 3(t-t') 
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Figure 2. Sampling of a chromatogram in the time and frequency 
domain (schematic). I, a continuous chromatogram and a delta function; 
II, sample of the chromatogram resulting from the multiplication of 
the continuous chromatogram by the delta function; III, chromatogram 
in reduced time; IV, Fourier transform of III; V, constant interval 
samples of III; VI, Fourier transform of sampled chromatogram, the 
dashed lines show the separate replicas of IV; VII, the Fourier 
transform of a sampled chromatogram that has less intense high 
frequency components than VI_ so that no overlap occurs, X 5g (8) is 
the ie, sampled form of X ,(s) Pe Miabdes v4 is the inverse 
Fourier transform of the frequency sampled data X 5d oR AY Re aN 


VIII the dashed lines are the separate replicas of the original 


continuous-data function in the time domain. 


LAG epind ‘ans 


} Sed) sap ERO g Bap sing 
PUI Bo anes 


ik am ! e 

He oR eT wea bus , bain te “kt a 

; inte: bay ett , hi Lanes epamanciig 

hg : . 1 csi oil t Milade se Lome wa mp oC 
'Py . 30 riled Wate OV" peekee 
Wise Le oe ee, Fe tapee oly ie 

rts a a i) { ‘ u 

vi ny SRO rer « FTX “Ay aad wot 


| aeil er nytseaast femoris pacensd 6 3 


\ 


imviy GR SRI We * Shp oc .. 


LEY, 45 i go Deities: 


S67 h099%, BISeRASe ONS of: 


iisneh stedosfl ani 1Garo 


1 : ( : " 
- ery ' 


! 


26 


where 
he) = oC) JAG 


and 2 is reduced time equal to t/At . The result of sampling the 
data, & (2) , is obtained by multiplication of two functions: the 
original data function and another consisting of a series of delta 
functions spaced one unit of reduced time apart. Representing the 
series of delta functions (the sampling function) by the symbol Il(z) 


in Equation 24 gives 


© (2) Str (27 ues) (25) 
where 
I(a) = v 6 (z-n) (26) 
In the frequency domain, 
X fs) = X(s) * Ils) (27) 


where s is reduced frequency equal to vAt , and X(s) is the Fourier 


transform of the continuous data (Figure 2, IV). 


Figure 2, VI shows the effect of convolution by [[[(s) in 
the frequency domain; replicas of the Fourier transform of the original 
data function, X(s) (Figure 2, IV), are formed at intervals of one 
unit of frequency. These replicas will overlap (Figure 2, VI) and 
interfere unless X(s) is equal to zero for all frequencies greater 


than or equal to 1/2. If there is no overlap, 4%_(s) will be 
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identical to X(s) over the region where |s| issless *thian 84/72 
and the inverse Fourier transform of the central replicate will 
regenerate the original continuous curve. Therefore, provided the 


sampling interval, At , is small enough so that 
X(s) = 0 when |s| > 1/2 (28) 


the samples will contain all the information in the data. In absolute 
units the critical sampling frequency 1/(2At) Hz is called the 


Nyquist frequency. 


The Fourier transform of a function sampled in the time 
domain is continuous (Equation 27), whereas the Fourier transform 
algorithm used in a digital computer is discrete; that is, the Fourier 
transform is calculated at frequencies that are multiples of 1/V 
where W is the number of samples taken. (WAt is equal to T , the 
interval over which the data are sampled in the time domain.) When 
the fast Fourier transform algorithm (14) is used, the number of 
samples must usually be a power of 2 , though the number of samples 
available may not be an exact power of 2. If M is a power of 2 
and the number of samples available, W , is less than M , then 
the last M-N samples can be given the value zero, and the Fourier 
transform will be calculated at frequencies that are multiples of 1/M 
rather than 1/W¥. The case of NW equal to M is discussed below, 


although the theory for the case of W less than M is different and 


more complicated, the practical consequence of sampling the Fourier 


MPeeformn 4s the same as that for N equal to, M . 4A discrete Fourier 


transform of a function is equivalent to taking constant-interval 
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28 
samples, in the frequency domain, of the continuous Fourier transform 
of the function. The sampling operation in the frequency domain can 
be represented symbolically as multiplication by a sampling function, 


N Il(Ns) (Figure 2, VII): 


X gS? = X fs) N Il(Ns) (29) 


The time-domain version of frequency-sampled data is the inverse 


Fourier transform of Equation 29: 


X eal? = x (a) * TI (2/N) (30) 


When the sampling theory used for sampling in the time domain is 
applied to sampling in the frequency domain, samples of X,(s) at 
intervals of 1/N are equivalent to the continuous function, Xs), 


provided the section of base line taken is large enough so that 


e(2) = 0 when l2| > W/2 (31) 


Figure 2, VIII shows the overlap error resulting from an 
attempt to regenerate a data function when WN (or fT ) is too small. 
Since a Gaussian function approaches zero only at infinity in the time 
domain, taking the discrete Fourier transform of a Gaussian peak will 
introduce an error the magnitude of which will increase as the width 
of the peak (9) increases relative to the period of time over which 


the data are taken (7) . From a table of ordinate values of a 


Gaussian function, the maximum error relative to the peak maximum is 


less than to when YZ is greater than 9.80 . Also, since the 


Fourier transform of a Gaussian function is itself a Gaussian function 
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29 
[with a standard deviation equal to 1/(21mc)] and is never equal to 
zero for any frequency, taking samples of a Gaussian peak in the time 
domain will introduce an error the magnitude of which will increase as 
the peak width decreases relative to the sampling interval (At) . The 
maximum relative error is less than Lie” when the peak width is 
greater than 1.7A¢ . The combination of these time- and frequency- 
domain restrictions means that at least seventeen (1.7 x 9.8) samples 
of constant interval are required to specify a Gaussian peak with a 
relative accuracy better than hor” i! The wnequirements of AG and dZ 
for a peak in the presence of base-line noise are discussed below in 


the section on estimation theory. 


The condition that a base line be zero at the ends of the 
section of base line over which samples are taken (Equation 31) is 
seldom met in practice because the variable zero control of the detector 
1] set so that the output signal will never be negative or zero 
throughout the chromatographic run. Also, trends in the base line, 
or noise with a seer eee than 7/2, will cause at least one end 
of the sampled section to be significantly different from zero. 
Although extra parameters for the mean and the trend in a base line 
could be included in the peak feats ft is easier to filter thesdata to 
remove these components before Fourier transforming. Including a 
parameter for the mean is equivalent to filtering the eet with a 
function that has zero area, but a filter that removes any component 


other than the mean may destroy some of the information in the data 


about the peak being measured. We have found that the noise-whitening 


(h(t) in 


filter used in the estimation of power-density spectra 
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Equation 20] is suitable for this purpose. 


Interval-Area Sampling 


Interval-area samples (13) are defined by 


abs 
eats ele) ari wre) (32) 
Equation 32 can be written 
© (2) =e ay ee) Tite) (33) 


where integration is represented as convolution by a rectangle function, 


iets) 


and 


, where 


(2) c=" when |2| Shy 2 


Il (2) 


H} 


O when, jal een/2 (34) 


In practice a small dead time occurs between interval-area samples while 


the integrator is transferring the data and being reset to zero in 


preparation for the next sample. If d 


is the relative dead time, the 


2t+d/2 


rectangle function in Equation 33 should be replaced by Il igre? , but 


if the dead time is small compared to the sampling interval, this small 


correction can be ignored. 


In the frequency domain, the interval-area sampled data become 


el Ke) Sint) * Ils) (35) 
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Comparing Equation 35 with the equation for instantaneous sampled data, 
Equation 27, interval-area sampling multiplies the Fourier transform of 
the original data by the Fourier transform of a rectangle function 
(Figure 3, I). This gives interval-area sampling an advantage over 
instantaneous sampling in that noise at high frequencies (Ss greater 
than 1/2) is attenuated and is less likely to cause overlap errors. 
The data in the useful frequency range (s less than 1/2), however, 
will be distorted. The determinate error in the estimates of peak 
parameters caused by this distortion can be calculated from the 
definition of cumulants (Equation 8). The resulting corrections are 
the well known Sheppard corrections for grouping (22); the corrections 
for peak area, retention time, and odd-order cumulants are zero, but 
the corrections for peak width and even-order cumulants can be 
Significant. These corrections could be applied to the parameter 
estimates or included in the peak model, with the understanding that 
neglecting high-order uncorrected cumulants introduces the prior 
information that these cumulants have the values of the Sheppard 
corrections. Sheppard's corrections are exact for a peak function, but 
only approximate when a peak is contaminated by noise. It is better, 
therefore, to correct the data rather than the parameter estimates even 
though more computation is required. Correcting the data, however, 


introduces a bias into the estimate of the power-density spectrum; more 


is said about this below. 


After filtering the Fourier transform of interval-area sampled 


data is 
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Figure 3. Digitization effects. I, Fourier transform of rectangle 
function;  I1, effect of jitter errors (schematic); III, power- 
density spectrum of jitter noise; IV, power-density:-spectrum of 

noise caused by integration dead time fluctuations; V, continuous 
data (schematic, with arrows representing instantaneous samples); 
weeeresult of the quantization of the data in V (schematic); 

VII, probability-density function of continuous data in V_ (schematic) ; 
VIII, probability-density function of quantized data in VI 
(schematic). The lined area in II is the interval area sample 
resulting from jitter. With no jitter, the integration would extend 
me thesdashed liness The *area of. the lined part of . VII. 1s equal to 
the area of the delta function represented by an arrow from the origin 
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X'(s) = X(s) H(s) outa). (36) 


where H(s) is the transfer function of the filter and where the 
sampling-interval and the sampling-time conditions (Equations 28 and 
31) are assumed to have been met. The corrected Fourier transform of 


the data required for the estimation of the peak parameters is then 


(Ss) Xi (s) [H(s) soni) (37) 


For the calculation of the integral in (2) Equation 15 we define 
P,A€ 
~ 


the Fourier transform of the data at zero frequency, X(0), to be zero 
because the reciprocal of the power density spectrum, one eS 


defined to be zero. 


Jitter 

Variation in the sampling interval is called jitter; as shown 
below, jitter adds noise to interval-area sampled data that increases 
with signal amplitude. Although we do not consider jitter noise in 
the noise model (because it is nonstationary), jitter always occurs to 
Some extent in practice. We consider jitter noise here with a view to 
measuring it and thereby determining the limitations of the noise 
model. The results below are a simple extension of the theory 
presented in (23) to interval-area sampling. Figure 3, II illustrates 
that the jitter at the beginning of integration, mo (e-l72) ey se coe 
difference between the time when integration actually begins and the 
& = 1/2, “Phe effectror 


time when it was expected to have begun, 


jitter is then as shown in the equation 
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BU5+E (2445) 
z..(a) = tf a(t) at) Il(z) (38) 


metaylor series expansion of w. (8) taken to the first power of 
the jitter gives 
teeta) ays fa) + [x(2+1/2) &(2+1/2) 
(39) 
setac- 172) atlerl/2) lena] 
Equation 39 is an adequate approximation of Equation 38 when the 
standard deviation of the jitter is small compared to the sampling 
anterval, andiwhen d2«(2)/d 2 is small compared to x(2) . The 
Slope condition holds when the Nyquist sampling condition is met and 
there is little signal density at frequencies near the Nyquist 
frequency. Jitter as defined above can be considered a continuous 
random process sampled at the same time as the data. It is assumed 
that the jitter has a normal p.d.f. with zero mean and standard 
deviation oe , and has a white power-density spectrum extending beyond 
the Nyquist sampling frequency so that the autocovariance function of 
the jitter is effectively a delta function. It is assumed, also, that 
Pieejyiiter, €(2) , fie whcorrebated wath the signal, <2(2)) 4@sif"the 
signal and the jitter are correlated (by mechanical vibration, for 
example) a more elaborate model than that used here would be required. 


The autocovariance function of jittered interval-area samples is then 


-—— yp, 
Rape (AP a (Oy se ) 0. keote) 
j — xe x J 
(40) 
-~§ (g+1)-6 (2-1) ] 


and the power-density spectrum is 
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owe? ~ Oy (8! + 2(W +a") 05 [l-cos (278) ] (41) 


The mean of the data, X , is included in Equation 41 because in the 
calculation of the noise power, Wo in Equation 22, the power- 
density at zero frequency is given a value that is reasonable for 
noise with a zero mean; the actual value of C,. (0) is equal to the 
Square of the mean of the data. Equation 41 shows that jitter adds 
a noise component to interval-area-sampled data that has an intensity 


proportional to the square of the mean of the data. 


Another noise source associated with sampling is caused by 
the variation in the integrator dead time referred to previously. 
Assuming that the dead time has a normal p.d.f. with a standard 
deviation, O53 , that is much less than the mean dead time, and that 
the variation in dead time has a white power-density spectrum 
uncorrelated with the signal, dead-time noise adds to the base-line 
noise a component with a white power-density spectrum equal to 


Wt) 0% . Figures 3, III and IV compare the power-density 


Spectra of jitter and variable dead time noise. 


Quantization 


The second step of digitization is quantization. As shown in 
Figure 3, V to VIII, quantization can give a data sample a random 
error and a constant error both of which may be as large as one half 


the quantization interval. Below the theories presented in (24) 


and (25) are summarized and interpreted. Quantization of a random 


variable is equivalent to taking interval-area samples of the p.d.f. 
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of the noise. In a normalized scale where the data, x » are taken 
in units of the quantization interval, the Dsdsfs of quantized 
noise is 

Pt) 2) p(x) * [[(~-1/2) [J (a-1/2) (42) 

q 

where p(x) Toecve. p.d.i,, OL the data at a specifica time, A 
considerable simplification in the equations given below results when 
p(x) in Equation 42 is shifted downward one half unit. Although 
shifting adds a constant error to the mean of the data, this error can 
be ignored since the mean is arbitrary. The joint p.d.f. of «x at 


i 


time t , and x. at time §(¢+r) is then 


20, BAe = Py (2i505) * W(x) * Ife) ([(x,) Ills.) (43) 


PaamtiesPourier transtorm of this joint p.d.f. (along the & axis; 


not the usual transform along the time axis) is 


sin(mx,) sin(mX,) 


= a7 Sea a ix? (44) 
Pe (Xp Xo) = P(X, Xo) ™X ™, * TxX,) * MX, 


The autocovariance function of the quantized noise can be obtained 
from Equation 44 with the formula 

eZ 

oP. (X4 1Xq) 


1 q eee (45) 
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When the data have a normal p.d.f., the Fourier transform of 


the joint p.d.f. of two data points is 
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4 ee 42 2 
PX 1X) = exp{ 9 IX #20 (TIX X5*X 51) , (46) 
where 


Ww = (20) “Wt ,/2 (47) 


In terms of the autocovariance function of the base-line noise 


(Equation 14), the variance, Ae , in Equation 46 is 
Wo = e,, (0) (48) 
ena thesautocorrelation function,,ip(t) » in Equation 46. is 


o(t) = 2, (t)/e,, (0) (49) 


From Equations 44, 45, and 46 the autocovariance function of 


quantized noise is 


° 2 
Corte = h Ott) + Ao * ) a exp (-wn ) 
x Cc = 
(50) 


- aan ) ‘ {22 exp{-w[n?+29 (x nim” } 


where the prime on } means omission of the value at nm=0. 
n 


Equation 50 shows that quantization adds noise to the data, 
and the second and third terms of the equation express the autocovari- 
ance function of the quantization noise. Using the symbol eae for 


the autocovariance function of the quantization noise and using only 


the most significant terms in Equation 50 gives, for the indicated 
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values of the autocorrelation function of the base-line NOLSS fe OAT, 


eZ (AW 41/0 ye" tere 
BURA ras UL (ee ar) pee Ot) y2 
oe aye” J) ert yao (51) 
ee eeewnenc) 2222, 0/2) saeae.t meee 
Ena oh (ai -1/n*)e” Ane Rey eRe 


Since the derivative of ade) wath) respect, to “pi t)seus 
cee ee One aawene, pit), as,equal.to.+ J... the, constants.(+ 2/712) 
meeeguation 51,, occur only when, p(t), is identically equal to, +1,. 
Pmmostapractical situations o0(1t) is identically. equal. to..+1 when 
ime tag is zero and is seldom equal to +1 when-the lag is any other 
value. Therefore, un? has a component at zero lag that is in 
effect a delta function. The power-density spectrum of quantization 
noise, consequently, has a white component that extends to infinite 
frequency, but when the data are sampled, only the portion between the 
negative and positive Nyquist frequencies is observed. (There is no 
overlap error because the autocovariance function, being a delta 
function, is not affected by sampling.) The impulsive nature of ¢_ (T) 
arises from the assumption that the quantization interval is exactly 
pene, jitter in the quantization interval no doubt occurs to some 
extent in practice and would broaden the impulsive component. The 
resulting power-density spectrum would approach zero at high frequencies, 


but unless the sampling interval were extremely small or the jitter in 


the guantization interval severe, the power-density spectrum of the 
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impulsive component should be nearly constant in the observed frequency 


range. 


Apart from the constants (+ 1/12) 4 the. ters in Equation 51 
do not vary rapidly with small changes in the value of 0 (tT); hence 
the equation can be used to sketch in the approximate autocovariance 
function of the quantization noise when the autocovariance function of 
the base-line noise is known. When the standard deviation of the base- 
line noise is small enough, the terms other than the constant 1/12 
can be neglected. (If the standard deviation is one half, in units of 


o 


the quantization interval, the largest of these terms is about Lo 
Therefore, under most experimental conditions the effect of quantization 
is to add white noise to the base-line noise. In absolute units, 

: : , ‘ 2 2 k 
the power of quantization noise is 1/12 q (volts ), where g is the 
quantization interval in volts. When the p.d.f. of base-line noise 
Heenormmal, the p.d.f. of quantized base-line noise will not be normal 
because the p.d.f. of quantization noise is rectangular. When the 


standard deviation of the base-line noise is several times the 


quantization interval, however, the effect of the rectangular p.d.f. 


can be ignored. 


Quantization of interval-area samples occurs after the 
integration over the sampling interval, so the correction for the 


distortion caused by interval-area sampling, Equation 36, should not 


be applied to quantization noise. In the calculation of the integral 
in p (£) in Equation 15, however, the same error occurring in the 
He 


is . 
data, Ja (v) |? , and the power density spectrum, Civ) , results in 


the cancellation of the error. Therefore, applying the correction for 
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41 
interval-area sampling to quantization noise will not affect the 
parameter estimates, but will introduce a consistent error in the 
estimates of the parameter standard deviations (defined below), which 


depend on C,(v) but are independent of the data. 
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ESTIMATION 


Information from the Data 

The object of estimation is to extract a peak from a base 
line so that the peak can be reconstructed as it would appear were 
there no base-line noise. Finding the peak that caused the data 
reduces to estimating the parameters of the peak model when the model 
is given. But, because the estimates will be derived from random data, 
the estimates themselves will be random, and it will not be possible 
to state the exact values of the parameters that caused the peak ina 
particular chromatogram. Knowing the probability structure of the 
noise, however, one can give the probability that a certain range of 


true values caused the peak in the data. 


The parameters of a peak on a base line are related to the 


data through the equation 
a(tsa) = fl(tsa) + g(t) (52) 


where f(t3o) is the peak function, g(t) the base-line noise, and 
a the set of parameters. For the purposes of this study, Qo is 
% 


the peak area, ay is the retention time, a. is the peak width, and 


a, and higher are the standardized cumulants, K, and higher. When 


a peak is not present in the chromatogram, the Hod, tar Oline late, 


me, 2s given by Equation 15; therefore, when a peak is present, the 
~\ 


eos. Of xs given Q (the true values of the parameters), is 


) = p,.fz-f) (53) 
%~% 


by the use of Equation 52. Equation 53 holds when the peak model, 
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f (tia) , represents the peak completely, or at least accurately 
enough SO that any deviation of the peak model from the true peak is 


indistinguishable from noise. 


Prior Information 

The true values of the parameters, ® may vary from 
experiment to experiment because of indeterminate variations in the 
experimental conditions (non-base-line noise) that may affect the 
parameters without having an effect on the base line. The knowledge 
or the lack of knowledge about the parameters for a particular 
experiment is best described by a p.d.f., PB, () , called the prior 

ay 

p.d.f. of the parameters (26) . Assuming that the experimental 
conditions that affect the true values of the parameters do not vary 
during the elution of the peak, one can speak of the "true" values for 
a particular experiment. (It is assumed, also, that the true values 
Of the parameters are practically independent of the base line.) The 
prior p.d.f. states the degree of belief in the hypothesis (26) that 
previously obtained values of the parameters still -app ly ‘torte 
experiment under study. In other words it expresses, in concrete terms, 


the state of the information about the true values of the parameters 


before the experiment is begun. The amount of prior information can 


have an important effect on the accuracy and precision of the parameter 


estimates; precise prior information is essential for the estimation or 


detection of small peaks that are barely distinguishable form the base- 


line noise. 


Often in an analytical problem the experimenter has some 


information, however vague, about the parameters before, the exper meat 
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is begun. In gas chromatography little is usually known about the area 
of a peak except that the area is finite and positive (or partly 
negative in anomalous cases), but the retention time and peak width of 
a previously identified compound are often known precisely. This 
prior information could be either obtained directly from the accumulated 
results of previous experiments on the same apparatus or predicted 
theoretically from published data and the experimental conditions. The 
prior p.d.f. of the parameters obtained by either of these methods 
is likely to be normal or nearly normal when the number of previous 
experiments is large by the central-limit theorem (11) . A normal 
prior p.d.f. is completely characterized by its mean, aa , and its 
variance, Ae PeaNGrinal tty OF the prior p.d.f. das unnecessary, gouc 
desirable because it aids the numerical determination of the most 
probable values of the parameters. If prior information is not 
available for a particular parameter, it is reasonable to assume that 
all values of the parameter within a physically possible range are 
equally probable; then the prior p.d.f. is uniform and can be 
represented by a rectangle function. An exact statement of the limits 
of the range of the parameter is not required. For the estimation 
method described below, it is important to know only that the true 


value of the parameter is well within the permitted range. 


The prior information introduced by the choice of the peak 


model can be interpreted as giving the nonadjustable parameters certain 


Peene and zero variance. For example, in the, Charlier derivative 


series with a Gaussian basis function, the high-order cumulants have 


means equal to zero. 
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Bayesian Estimation 

Erom'the “pidst.- %of “the data given the parameters 
(Equation 52) and the prior p.d.f. of the paraneters) therto7dcry 
of the parameters given the data can be obtained from the following 


form of Bayes' theorem (27): 


Py (a) Pool IR) 
Pala le) - a (54) 
~% 
where 
P(e) = op. (w) Poi, Zig) ag (55) 
%—v ~% NON 


The multiple integration in Equation 55 is over the parameter space 

Q . Although not indicated in the notation, Equation 54 is conditional 
on the experiment being run in situations where the prior p.d.f. of 
the parameters and the previously obtained autocovariance function of 
the base-line noise are applicable. The p.d.f. of the parameters 
given the data defined by Equation 56 is called the posterior p.d.f. 
of the parameters (26). The posterior p.d.f. does not contain 
irrelevant information about the details of the base-line noise but 
does contain all the essential information in the data about the 
parameters. Being a multivariate function defined over an infinite 
domain, however, a posterior p.d.f. 
Since numbers are easier to deal with than functions, a few character- 


izing numbers will simplify the interpretation of a posterior p.d.f. 


Below we discuss the estimation of parameter values that best describe 


the position or location of a posterior p.d.f. 


is too cumbersome for most purposes. 


ertyotet ota ioe Ror dndie, ae 02 ‘heh hail fines « 


“gnay 


(aa) 


; susbnfa ~etoménay ailz weld et BE seems re: voontinaael 4 


Lanai tomas a he seen puoi sith ing asegolgih: 4 


She 


te noOLsore aome bsevtimnir piabtnards valk uae dy bane 
y Paptone ses - 20 


a ib. “* sola dgean ‘eas batten: ah: we ‘hethaaparet, yi, balstae, ch § 


‘spega exon rot 


“ axegoneseg, seit rhea oh e ie. ae . 5 a 
aS 2, a eae eros oe ult 20: eke sepa. a | 


yous awh viet spoke th ire Sse 


iass9es Seat eae | 
» ae gaion enttt—and aula shiganb tits ets Rapa 
ped won’ sls $f A pestegaehan tek dsQaem okt, Ke 
sectqar tt ad ape Hien debdomet stnisevizing ® poked” 


ee 7 im AG : 7 » i 
y " r Pee ' 
i fh yn 


Cay’ moncbit | ‘ wage 


sbi. ck: ak aria iiiiialedl otis win 


ee ved ae rt? kage: 


mee ed «TF; fee Isirsseqy # pane . 


“epadsiranily 8? de amet it, iw tae 0 vo kave ow eto 


8 Bo “yoauteed & Bo wodandsath sir atts GHamte Liv 
nh et ae dint et vane 


a 
| : : : 7 
wt “= y Be, +4 > : 

roe > — 
¥ , ; : 
7 oR 1s —_? 
: . 
= 
i i. | —. 


serartet valcacdong 6 je idindot wa vo. 


46 


Point Estimators 

A set of operations applied to data to obtain an estimate is 
called an estimator. The resulting set of estimates is represented 
by the symbol a to distinguish it from the set of true values, a 
Our goal is to find an estimator whose estimates contain all of the 
available information. Such an estimator we call "ideal". Since the 
posterior p.d.f. contains all the information about the parameters, 


the relative merit of an estimator will depend on how well the p.d.f. 


of the parameters given the estimates, p iq Re) , "matches" the 
ny) 
7 L] 
pectexrior p.d.f., Pale GIR) » Ih the two p.d.t." sare not tne 


seme, the estimates willl. contain less information than the data. 


When the moments of two p.d.f.'s are the same, then the 
p.d.f.'s themselves are the same; therefore, moments should be a 
valid basis for the comparison of a posterior p.d.f. with an 
estimator p.d.f. The first two moments are usually sufficient to 
characterize an estimator. It is conceivable, however, that if some 
estimator A were poorer than estimator B on the basis of the first 
two moments, estimator A might "match" the posterior p.d.f. better 
than estimator B on the basis of higher moments. In practical 
applications characteristics of an estimator other than moments may be 
important. Some examples of practical criteria are economics, 
Mathematical simplicity, calculation speed, and insensitivity to changes 
in the probability structure of the system (in particular, nonstationar- 


ity in the autocovariance function of base-line noise). These criteria 


are not considered here on the grounds that once an ideal estimator is 


obtained, a simplification that results in a more practical estimator 
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can be justified by comparing the moments of the practical estimator 
with those of the ideal estimator. Furthermore, the mathematical 
structure of an ideal estimator may contain clues to practical 
simplifications; working downward in complexity is easier than 


building a better estimator from a poor one. 


The difference between the first moments or means of an 
estimator p.d.f. anda posterior p.d.f. may be called the bias. 


The bias for parameter y is 
|x) da (56) 


By use of the expectation operator, Ey , to represent averaging over 
Peetehat as, multiplying by the p.d.f. of variable” y “and” integration 
with respect to y ), Equation 56 can be written: 


Bah oy = 
x. 


toe = Bansal) (57) 
J 214 | 


a 
J uy 
The second moments about the means or the variances are 


{[a.-E (a 5 (58) 


“a; Fale SF “ale 


“Gi, 
os) 
and 


2 
)ee Alo.-E,,  (@,)1 } (59) 
‘a; ® ~ Male “Fala 


An ideal estimator will have zero bias, and the variance of the estima- 


tor p.d.f. will be equal to the variance of the posterior p.d.f. 
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Although an evaluation of the bias and the variances defined 
above may be difficult, it is possible to approximate them or at 
least make inferences about them by considering averages of the 
estimates over the data and the parameters. Averaging over the data 
for a given fixed set of values of the parameters, a conditional 


bias can be defined by 


ch neues (60) 


AV) 


cil = "x 


The conditional bias is an average over the data of the difference 
between two parameter values, while the bias of an estimator defined 
in Equation 56 is an average over the parameters of the difference 
between two p.d.f.'s. In some situations using an estimator that 

has zero conditional bias may be desirable, but as shown in an example 
below, an ideal estimator may have a non-zero conditional bias. The 


overall bias of an estimator, averaged over both the data and the 


parameters, is 


be =f 5 61 
B. E ee eee (61) 


where the joint p.d.f. of the data and the parameters Fequired an 


the expectation operation, E, . , is 
KR 
= 62) 
(7,0) =p (a) p. (|) ( 
CE a ae 


The overall bias also is different from the bias of an estimator 


defined in Equation 56. If one is zero, the other may be nonzero, 


but under certain conditions discussed below, both will be zero. The 


parameter a. is 


conditional variance of an estimator for j 
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2 
V.(a) =E {[@-F (257 63 
mee alan ge a 
and the overall variance is 
Vos Can We ic) (64) ° 
J Fit J iM 


The overall variance can be obtained experimentally by running a large 
number of experiments; usually, it is the only characteristic of an 
extimator that can be directly observed in practice. The other 
characteristics must be inferred from a p.d.f. that may not be 
directly observable when the p.d.f. is a function of the unknown 
true values of the parameters. When an estimator is unbiased 

(overall bias equal to zero), the overall variance is the variance of 
the estimates about the mean of the prior p.d.f. oe , in the case 
mea normal prior p.d.f. .and the midpoint of the range in the case 
Breasunizorm p.d.f£.). Since the overall variance is strongly influ- 
enced by the randomness of the parameters, it is more characteristic 
of the entire experimental system than of the estimator. Because it 
is less influenced by the prior p.d.f., the mean-square error, defined 
as the average difference between the estimated and the true values of 
the parameters, characterizes the estimator better than the overall 
variance. The conditional mean-square error is 


2 
= aa. (65) 
5 .(g) = Baal ( ae 


and the overall mean-square error 1S 


2 
a ae (66) 
S. = E— [(a, o2) ] 
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The following equations give the relations between the mean-Square 


errors and the biases and variances: 


\R 2 
o (a) = tat + ey (67) 
Se =.= VY. + °2E -B 
j j oe 9 (LB B(o)]} (68) 
J 
where YV is the variance of the prior p.d.f. The conditional 


Pj 


mean-square error is always at least as large as the conditional 
variance. The overall mean-square error may be larger or smaller 
than the overall variance; if the conditional bias is zero, it will 


be smaller. 


An important result concerning the overall mean-square error 
Gan be derived from the posterior p.d.f. This result, called the 
Cramér-Rao inequality (28), states that the smallest possible value 
of the overall mean-square error of any estimator is greater than or 


equal to the diagonal of the dispersion matrix: 


Oe Se Thi 6 (69) 
thse | 


where the dispersion matrix is the inverse of the information matrix, 
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The elements of the information matrix are derivatives of the posterior 


p.d.f. averaged over the data and the parameters: 
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] C7a} 


The Cramér-Rao inequality holds provided the derivatives in Equation 71 
exist and are absolutely integrable with respect to £ and a . In 
the noise and peak models used here these conditions will not be 
Satisfied when the prior p.d.f. of a parameter, a, for example, is 
uniform because the derivative with respect to oy, Wil notVexisbuat 
the discontinuous end points of the p.d.f. The mean square error 

of oe might then be less than soa {BUCY ifethes posterior  prd.£, 
is concentrated well within the ends of the uniform p.d.f., then the 
Cramér-Rao inequality will hold approximately. An estimator that 
Satisfies the equality in Equation 69 is said to be efficient, 
efficiency, = , (overall mean square efficiency) being defined as 
the ratio of the minimum mean-square error to the actual mean-square 


CELOT: 


Steck (72) 
5 5g 3 


An efficient estimator exists only if the information matrix, defined 


by Equation 71, is independent of the values of the data and the 


parameters; this is true only when the posterior p.di.f. , is; normal 
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MP Estimator 


From Equation 57, the mean Of the: ~p.d.1.. of anwidedgt 


estimator and the mean of the posterior p-d.f. should have the same 


value. Therefore, the mean of the posterior p.d.f. would be an 
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ideal estimator in terms of bias but the mean is often difficult to 
calculate. Another estimator, also based on the posterior Deen We 
more eaSily calculated than the mean, consists of taking the most 
probable values of the parameters to be the parameter estimates (MP 
estimator). If an efficient estimator exists, the MP estimator will 
be efficient and unbiased, and the p.d.f. of the MP estimator will 
Ppemenessame as the posterior p.d.f. Consequently, if the posterior 
p.d.f. is normal, then the MP estimator will be ideal in the sense 
that the estimates along with the information or dispersion matrix will 


contain everything that can be said about the parameters. 


Example 

To illustrate the use of an MP estimator, consider the 
Provlem of estimating the mean, a , of the retention time of a gas- 
chromatographic peak from m measurements (runs) Ly oreer ik, s Runs 
on previous days showed that the true mean has a normal p.d.f. with 
Mean a and variance oe , and that the individual measurements are 
independent and have a normal p.d.f. with variance Go athe 
posterior p.d.f. of the true mean, a , given the data, x , is 
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(alae) = (2nV") S/20n (73) 
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The MP estimator, a, is found by equating the first derivative with 
respect to a of the right side of Equation 73 to zero, replacing a 


metre ad , and solving for a 


Ve2x+aC/n 
V +C/n 


The MP estimator of the mean, the weighted average of the prior mean 
and the data mean, is not the same as the commonly used estimator, x ; 
nevertheless, a and x become approximately equal when either the 
number of measurements or the prior variance becomes large. The 


Setimator p.d.f. of a, given a, is 


Pala (ela) = (an!) */? exp[-(a-a)*/2V"] (75) 


i i a c all values of a the 
Since Pala (ela) is equal to Pal (Ol®) or ' 
MP estimator is ideal; also, since the posterior p.d.f. is normal, 


the MP estimator is efficient and unbiased. 


The different quantities defined above that are used to 
characterize an estimator are listed below for the MP estimator of 


the mean retention time. 
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iL 


— 
role 
I] 
S 


TeC7 (nV) (78) 
p 
1+ (a =a) “C/ (nv) 
p 
£ 
ie acy an (80) 
p 
a wary: ae C/ii 
es = yeh es 07 ee (81) 
Pp 
V0 (®) = V6) = J' (82) 


The overall bias B and the estimator bias B..(&) are zero, whereas 
the conditional bias B(a) is nonzero except when the true mean 
happens to equal the prior mean (Equation 76). But, when m or as 
are large, the conditional bias approaches zero. A nonzero conditional 
bias would be undesirable when the estimates of the retention times of 
two peaks on the same chromatogram are used to calculate a retention 
matcioO. Or a retention index. Here, it would be better to use the MP 
estimator of the final result, utilizing prior information of the final 
result itself, rather than to combine the estimates of two retention 
times. The estimator variance, V (a) , is less than either pines prior 
variance, , or the variance of the average of the data, x , which 
is equal to C/n . The overall variance, JY , is greater than the sum 


of the prior variance and the variance from the data (C/n) when C/n 


is less than V ; in other cases JV is less than Vs (Equation 80). 
p 
If the prior variance was assumed to be infinite when in fact it was 


finite, that is, if the prior information was ignored in the estimation, 


the overall variance would be equal to the sum of the prior variance 
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55 
and the variance from the data, i + C/n . In the statement of the 
estimation problem the true mean was assumed to be random; in some 
Situations the true mean may be constant and the prior information 
may be random. For example, suppose the prior mean, a , is the 
average of m measurements and the prior variance, ke , iS equal to 
C/m . The MP estimator would in effect be the average: Of) 7s plus 
m measurements, the m measurements being considered as prior 
information. The overall variance (Equation 80) correctly gives the 
average variance of aq about oe . But, because the prior mean could 
be far from the true mean in a single sequence of m runs, the overall 
variance is not a useful result in this case. The overall mean- 


Square error gives the more realistic variance of the estimate about 


the true value of the mean. 


Estimator for Nonlinear Peak Model 


In the noise model described here the parameters must be 
linear in the peak function for the MP estimator to be efficient. 


Edgeworth's form of the GCA series (Equation 10) is nonlinear 


(except for a the peak area), and therefore the MP estimator will 


O 7 
be inefficient and generally biased. Nevertheless, the loss in 


efficiency and the bias are not usually serious. The primary hypothesis 


of this thesis is that the MP estimates along with the dispersion 


matrix or the information matrix, though not ideal characterizations 


Sieasposterior..p.d.f., are approximately ideal. Even when the MP 


estimator is poor in terms of variance and bias, it may be valuable 


in its own right in giving the estimates that are the most probable in 


the light of the data available. 


i - 

whe th tacit lows it , NG + Y- seb off dem aimee. * 
state al gabigs |a& of nana Ger igh. su a7 orld nba - ya 

- on =, PAD “inte Ps tH tad Bf "4 ar | Vii uyom ie in ; 

1 au : x 
Wf encod .siqeetie TOT: : 


we Sab oi % ties vealed sad hms. amore 


y = {> AE 7 tn ; gs alex Jeni ties GA: ; 


rie sah mn r¢ si b iis os nel ea q 
' % Yo anhMighas sien — “ geen «oes geet wit 


+s } i 'Aksraha coffers £ son ee z 
' ‘ec¥ Sc 5 <t > Jideett. eaves 


wel 
nase of ie teks Sr 


in 


Ne * 
LarOw 5 ad LCgtik. 5 thay or 


eee 


' 
SUM jee f ne 682 ere Soetiess Luba Setar a at 
\ 7 


cas ooy ,oeemisgear?® git “ea} ee) nofecist Seed eae 


« 


: . 4 P 
: ft 
aishhinon =) "(Ci noalpar emgree “Att oi 34 went = 
iy wttinites Gu. cats gopteweds rs. .(hotkh veeg of . s” wot 348 
- e ‘ 
_ AE emet ‘at ,weateteyeveG .feanit yiicussep See Je 


cteatingy! Seattey Sty .Randdde Gidelin: 2054 <r ccld =e, Bee 


aT 


gotetegath SA MER Bela Cee es «<A tS feds el ated ; 


‘ 
eanisezr | veh tals abi jag iqyats ,4is3io robiaswin wi 30 c.3 33; 
in 


7 


. ~ Fi 
 ditinkfiw nowd =. feels gietanieecigge crs . 2.8.5 
7 _ i 


sitnadat od “eet 21 . ana Tere Gahdtee t> =nva? o: 2edg aS 


. ab tidaigeg teen ich Axe 2) Ravendsey ol: yorvee at eégis. 
| - Wat 


56 

The Edgeworth series is a transcendental function of the peak 
parameters. Hence the MP estimators for the individual parameters 
cannot be derived in the form of an explicit algebraic function as 
can the MP ee eietor of the mean in the example considered above 
(Equation 74). The estimator of each parameter depends* on the values 
of all the other parameters, and all the estimates must be determined 
Simultaneously. If the prior information about a particular parameter 
is strong and the information from the data adds little to the total 
information (see below), simply assigning the value of the prior mean 
to the estimate of that parameter may be adequate without calculating 
it, that is, the parameter may be made nonadjustable. When only the 
first few parameters of a series peak-shape model are estimated, it is 
assumed that the data will contribute little to the information about 
the high order parameters; these parameters can thus be given a constant 
value and made nonadjustable. When the peak-shape function is well 
known, it may be necessary to estimate only the area, the retention 
time, and the width of the peak. Reference (15) gives an example 
in which the parameters of the peak-shape function are a set of samples 


of the peak on a previous chromatogram. 


Fisher's Measure of Information 


The posterior p.d.f. (Equation 54) is proportional to the 


meoduct. of. the prior .p.d.f. and the p.d.f. of the data, given the 


values of the parameters. Therefore, the information matrix defined 


by Equation 71 is the sum of two parts: 


“the dependence is mathematical rather than statistical, 


though the estimates may be statistically dependent also. 
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The prior-information matrix, fd , has the elements 


2 
dein P(e) 
] (84) 


if the vector of prior variances, te , 1S extended to be a matrix 
with covariances as well as variances, then the prior-information 
matrix will be equal to the inverse of the covariance matrix. When 
one or more of the parameters has a uniform prior p.d.f. and if the 
MP estimates are well within the end points of the uniform p.d.f.'s, 
the corresponding rows and columns of the prior-information matrix 


can be given the value zero. The data-information matrix, JI, , has 


va 
the elements 
Agee re! Po a B18) 
a See eer nS ee (85) 
ak x gun 3 a7, =A 
In the noise model used here, Equation 85 becomes 
S Te» 9 Eevee F' (s) 
7 sa f at os a ds (86) 
: S) a: 
dak -o “2 J k 


where the dagger +t means take the complex conjugate. 


The diagonal elements of an information matrix are quantitative 
measures of the average amount of information contained in the data 


(including the prior data) about the parameters, or from another view- 


point, they are the worth of the parameters in terms of the information 
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the parameters contain about a peak. Some parameters may contain 
information from the data already contained in other parameters. The 
off-diagonal elements of an information matrix express the degree of 
this redundancy between two parameters. Although there are other 
measures of information, this one, called Fisher's measure (30), 
arises naturally from considerations of the posterior p.d.f. and is 
intuitively related to the familiar usage of precision defined as the 
reciprocal of variance. Fisher's measure of information is not an 
absolute measure of precision or information, however, because the 
Magnitude of the matrix elements aoe depend on the units of the 
parameters. Therefore, Fisher's information measure may not be suitable 
for the comparison of the information content of the parameters of 
different peak models. The magnitudes of the diagonal elements of the 
information matrix are used here for the comparison of the information 


in the parameters of peaks with different parameter values. 


An essential requirement of a series model of a peak-shape 
function is that the information content of a parameter in the series 
must decrease with increasing order of the parameter. As more and more 
parameters are included in a model, they should eventually contain all 
the information about the peak*. Edgeworth's series satisfies this 


requirement: in white noise the ratio of the information in parameter 


2 
r to the preceeding parameter of order (r-1) is equal to 2(2r-1)/r 


when r is greater than 3 In general, the information matrix for 


any peak model is symmetric, element Pia is equal to element Aue : 
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25) 
and the determinant of the information matrix is greater than zero. The 
off-diagonal elements of the information matrix for Edgeworth's series, 
because of the symmetry with respect to frequency of the derivatives 
in the Be finition of the matrix elements (Equation 85), are zero for 
even-odd subscript pairs except for the elements containing derivatives 
Oeetne peak width, a... That is, element ae is equal to zero when 


2 


g +k is odd except when j or k is equal to 2. 


Predicted Standard Deviations 

The square root of the diagonal elements of a dispersion 
matrix (Pp , Equation 70) are estimates of the standard deviations of 
the parameters. These predicted standard deviations are random 
variables because, in the case of inefficiency, they depend on the data. 
The off-diagonal elements of D are covariances. They indicate the 
Magnitude of the statistical interdependence of errors between two 
parameters in the same way an autocovariance function characterizes 
noise. For example, if the noise were to cause parameter estimate 
@. to have a certain error, then parameter ay probably would have 
an error fn, times the error in a . The degree of interdependence 


‘ . ° ‘ * r) 
is more easily seen in the correlation matrix, ry , where 


1 
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(87) 
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The diagonal elements of a correlation matrix are equal to 


off-diagonal elements can be neither greater than “i? nor less* than 


-l. The Cramér-Rao inequality does not give specific information about 


the magnitude of individual covariances; the actual covariances may be 


greater or less than the corresponding element of the dispersion matrix. 
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A geometrical interpretation of the Cramér-Rao inequality in reference 
(28) uses the complete disperson matrix to make probability statements 


about the errors in the parameters taken as a group. 


Ideally the off-diagonal elements of the correlation matrix 
for a peak model are all zero. In this case the estimation of the 
high-order parameters of a series model could be omitted without 
Significantly affecting the estimates of low-order parameters. 
Neglecting a cumulant in Edgeworth's series, or giving it the value 
zero, when in fact it has the value XK!’ , introduces a determinate 
error equal to -K' . If lower order cumulants have a nonzero 
correlation with the neglected cumulant, the error can propagate 
downward to give the lower cumulant estimates a significant bias. 

On the other hand, if the estimator is inefficient, zero correlation 
does not necessarily rule out this kind of error. For Edgeworth's 
series the correlation between retention time (an odd-order parameter) 
and any even-order parameter is zero. Hence, neglecting even-order 
cumulants should not have a significant effect on the bias of the 


retention-time estimate. 


The other off-diagonal elements of the correlation matrix 
for the parameters of the Edgeworth series are generally nonzero, and 
their magnitudes vary greatly with the values of the parameters. The 


values of even-even and odd-odd subscript correlations are positive 


and may be an high as 0.99, but decrease as the. distance from the 


matrix diagonal increases. The values of the even-odd correlations are 


generally small or negative; sometimes they may be as low as -0.99. 
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Efficiency and Bias 
For an MP estimator to be efficient, the magnitudes of 

the elements of the information matrix must be constant over all values 
of the parameters. Therefore, as the values of the elements become 
more nearly constant over the range of parameter values in which the 
Posterior ipsdsf.ocis most«significant (such as + 3 standard 
deviations from the maximum), the MP estimator will become more 
efficient. As the number of independent samples used in the 
estimation approaches infinity, the MP estimator becomes efficient, 
unbiased, and normal (31). But the number of independent samples that 
contain information about a peak is limited by the finite width of the 
peak and by the autocovariance of the noise. Nevertheless, as the 
amount of information (the magnitude of the diagonal elements of i ) 
increases, the posterior p.d.f. becomes more concentrated about its 
maximum and the range over which the posterior p.d.f. is significant 
becomes smaller. Therefore, since the nonlinearity of the system that 
contributes to inefficiency does not increase with the values of the 
diagonal elements of a ; the posterior p.d.f. must become more 
nearly constant over the significant range, and the MP estimator must 
become more efficient as the information increases. This argument is 
often stated in terms of signal-to-noise ratios (32) , but these are 
equivalent to the diagonal elements of the information matrix, f. - 


In a situation where the MP estimator might be highly inefficient 


and highly biased because of insufficient information in the data and 


little prior information, to present the entire posterior p.d.f. would 


be better than to attempt to characterize it with a point estimate and 


an approximate variance. The amount of information about the parameters 
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can be increased experimentally by increasing the area of the peak 
relative to the noise power or by decreasing the similarity of the 
peak to the noise waveforms. For example, a chromatogram with a noise 
power-density spectrum concentrated at low frequencies would give 
more information about narrow peaks than broad peaks of the same area 
because narrow peaks have more information at high frequencies, where 


the noise intensity is low. 


Bias, Since it is defined in terms of means, can result from 
choosing the maximum rather than the mean of the posterior p.d.f. 
when the posterior p.d.f. is asymmetric about its maximum. Bias in 
an MP estimator caused by asymmetry in the posterior p.d.f. will 
decrease with increasing information for the same reason that 
efficiency increases with more information. Nonstationary base-line 
noise can cause bias in the estimates of parameters as well as. in the 
estimates of parameter standard deviations. Nonstationarity that 
results in changes in the magnitude of the power-density spectrum 
but not its mathematical form might be reduced by approximating the 
power-density spectrum with a simple function and estimating one or 
two parameters of the function along with the peak parameters. This 
method can be used only when the power-density spectrum does not 
change or changes only slightly during the elution of a peak. In the 
problem of estimating the mean of a set of measurements discussed 
above, the variance of a single measurement is assumed to be 


Constant and known. If the variance changed from day to day, it 
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could be estimated along with the mean. However, depending on the 
form of the prior p.d.f., the MP estimator of the variance may be 
inefficient (because the posterior p.d.f. is not normal) and may be 
biased (because the posterior p.d.f. is asymmetric about its 


maximum) . 


Sampling Interval and Sampling Time 

In practice an important cause of both bias and inefficiency 
is presenting the estimator with insufficient data either in the 
frequency domain (sampling over too small a range of frequency) or in 
the time domain (sampling over too short a period of time). Consider 
the loss of efficiency, and also the loss of information, in the 
frequency domain caused by sampling the Fourier transform of data 
between finite Nyquist frequencies. According to Equation 87 a 
diagonal element of a data information matrix, te. , iS an integral 

Jd 

ef a function of frequency. This function can be called the 
information function of parameter oe . The information functions are 
real, positive, and even functions of frequency, and they approach zero 
at high frequencies. The total information is the integral of a 
frequency function over infinite limits. Although using finite 
integration limits will cause loss of some information, a small loss, 
such as 1% of the total information, may be acceptable in practice. 


or some other acceptable information loss, integration 


famuts +S. .could be »found so that integration of the information 


Biieeion between ichese. limits would. attain, (100-1)% ,of, the total 


information. If the integration limits are evaluated for each parameter 


and the largest limit (s,) found, then the sampling interval At, 
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equal to 1/(28_) would be the largest that could be used for the 
parameters to contain at least (100-1)% of the information in the 
data. Outside the frequency range pice the data contain little 
information Ebout the parameters. Since data that contain no 
information cannot contribute to bias, the maximum sampling interval, 
At. , will be small enough to ensure that bias as well as information 
loss are within acceptable limits. Similarly, the time intervals, 

2 , within which the data contain an acceptable fraction of the total 
information about parameter oe can be found by use of the inverse 
Fourier transforms of the information functions. The inverse Fourier 
transform of an information function is a function of time and its 
integral between infinite limits is equal to the total information, 
jg t 3 However, finding the shortest acceptable time interval, i . 
Jd 
in the time domain, where information functions are oscillatory, will 
be more difficult than finding 3. in the frequency domain. The 
information functions in the time domain are symmetrical about the 
retention time of a peak. Therefore, data should be taken the same 


distance (2/2) on each side of the mean even when the peak has a 


pronounced tail. 


Even though the above methods for finding the maximum sampling 


interval and the minimum sampling time are exact and clearly defined, 


they may be difficult to evaluate. An approximate method for calculating 


the maximum sampling interval is discussed in the section on digitization. 


Points on a base line separated by more than one half the autocorrelation 


Width of the noise, T (Equation 22), are effectively independent; 
a 


therefore points more than one half the autocorrelation width plus four 
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peak widths (standard deviations) from the retention time of a peak 
should contain little information about the peak. Hence, a reasonable 


approximation of T. is (8o+2',). 
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CHAPTER 2 


EXPERIMENTAL 


Many computer-based methods for the reduction of chromato- 
graphic data are designed for the computer rather than for the 
chromatographic system. Most computer methods are designed to meet 
the limitations of data storage, time, and cost, or to handle data 
from several chromatographs simultaneously. These methods have an 
important role in the automation of routine analyses; nevertheless, 
if all the information in a chromatogram is to be used to obtain the 
most precise measurements, a method designed for the chromatogram is 
required. In research applications and in some routine applications 
the information contained in a chromatogram about a peak may be a more 
important consideration than the cost of extracting the information. 

The maximum postertor probability (MP) method, described in Chapter 1, 
is suited to such applications because it utilizes as much of the 
available information as possible for the estimation of peak parameters. 
In this chapter the MP method is applied to the estimation of peaks 
obtained from an isothermal gas-chromatographic system with a thermal 
conductivity detector. Isolated peaks are considered primarily, 


though in limited circumstances the method is useful for the detection 


of small peaks under a larger peak. 


In the MP method the most probable values of the parameters 


of a peak model are taken to be the estimates of the parameters of an 


experimentally observed peak. The probabilities of different values 


of the parameters are derived from a probability-density function 
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67 
(p-d.f£.) called the posterior p.d.f. of the parameters. Once the 
posterior p.d.f. has been obtained, the determination of the values 
of the parameter estimates is a purely numerical problem. The 
mathematical details of the method used by the authors are given in 


the Appendix. 


The posterior p.d.f. contains all the available information 
about the parameters of a peak. The MP estimates and the minimum 
mean-square errors calculated from the posterior p.d.f. together 
contain nearly all the available information in that they approximately 
specify the posterior p.d.f. The square root of the minimum mean- 
Square errors, here called the predicted standard devtattons, are 
approximate standard deviations of the parameters due to base-line 
noise alone. These and the tnformatton matrix, which provides a 


quantitative measure of information, are defined in the Appendix. 


The posterior p.d.f. is proportional to the product of two 
ped. f.'s (Equa onpA- 18 Appendix)! mOnel 1s thei) p: dat @oki the 
parameters given the data, which contains the information in the 
experimental observations about the parameters. The other is the 
prtor p.d.f. of the parameters, which contains information about the 
parameters available from previous experience. The p.d.f. of the 
parameters given the data reflects the effect of the randomness of the 
data on the apparent values of the parameters. The derivation of this 
p.d.f. requires clearly and carefully defined models for the randomness 
of the data, for the peak as a function of the parameters to be 
estimated, and for the relation between the peak model and the data 


(see Appendix I). These models may be called the a prior 
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information about the system because they are not based on direct 
experimental evidence. In contrast to the qa prtort information, 

that contained in the prior p.d.f. of the parameters may be obtained 
from direct observation. In this Study, prior information 2s 
unavailable, and therefore a prior p.d.f. is chosen that is the 
Beast informative,of all possible p.d.f.'s. A prtort information 
about a chromatographic system is obtained by a process of induction 
from other systems that the experimenter believes are similar to the 


system under study. 


Since a prtort assumptions may be incorrect for particular 
experiments, the results may not be valid descriptions of reality. 
Nevertheless, the results will be self-consistent with respect to the 
a prtort information, and independent checks to confirm the reality 
of the models being used can be made after the experiment has been 


carried out. 


For the purpose of obtaining all the information in a 
chromatogram about a peak, the mathematical form of the model for the 
peak is unimportant. The primary consideration is that the peak model 
must specify the actual peak as completely as possible. The peak 


model used here (Edgeworth's form of the Gram-Charlier type-A series) 


is based on the Gaussian function. The choice of the parameters used 


to indicate deviations of the peak shape from the Gaussian form also 


is unimportant. The shape parameters used here are standardized 


cumulants. (The word standardized is usually omitted.) Cumulants, 


though related to moments, are mathematically more simple and for this 


reason are preferred. The parameters are represented by the symbols 
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A for peak area, vo for retention time (mean), o for peak width 


(standard deviation), and K, Co) Kk for the cumulants. 


6 


In the model for base-line noise the average behavior of 
base-line noise is completely characterized by an autocovartance 
funetton, the Fourier transform of which is a power-density spectrun. 
It is assumed that the noise has a normal p.d.f. It is also assumed 
that the power-density spectrum is independent of time (stattonary) 
so that the spectrum can be estimated from a section of base line 
without a peak. Although the estimate of the power-density spectrum 
may be considered prior information about the chromatographic system, 
the assumption that the estimate is stationary and is an accurate 
description of the noise during the elution of a peak is a prtort 
information. Only noise that adds to the chromatographic signal is 
considered in the estimation of peak parameters. Multiplicative 
noise or noise with a variance that increases with increasing signal 
smplitude was studied and found to be negligible in the chromatographic 


system examined. 


In practice, base-line noise is not the only cause of 
variation in the values of estimated parameters. Fluctuations in the 
carrier-gas flow rate and the column temperature cause variations that 
may be larger than those caused by base-line noise. For the purpose 
of testing the MP estimator a chromatographic system is required in 


which only base-line noise affects the estimates of the parameters. 


Moreover, the true values of the parameters must be known to detect bias. 


A system satisfying these requirements to a degree that can show the 


influence of the estimator itself on highly precise results may be 
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Pier tculteto. constructs Such a system, however, can be simulated by 
placing peaks generated by a computer with known values of parameters 
onto sections of actual base-line noise. Although the simulated 
chromatograms can be used to test the MP estimator under the 

@ prtort assumptions, they cannot be used to test the applicability of 
the a prtort assumptions to reality; for this, actual peaks must be 


analyzed. 


In summary, the application of the MP method to a specific 
chromatographic system consists of the following steps. 
1. Characterize the random part of the data. Choose a model, determine 
the values of the parameters of the model (here we determine the power- 
density spectrum), and then test the applicability of the model (for 
Stationarity, normality, etc.). 
2. State the mathematical model of the peak(s). 
3. State the prior information about the parameters of the peak in the 
form'of a p.d.f. 
4. Write computer programs and test them with computer generated peaks. 
5. Analyze actual chromatograms. 
6. Study the results ‘to check the applicability of the models and the 
adequacy of the digitization of the data in steps 1 and 5 (quantization 


level, sampling interval, and sampling time). 


In the material that follows, presented in three separate 
parts, base-line noise is studied through power-density spectra, the 


MP estimator itself is tested through the use of computer-generated 


peaks on actual base-line noise, and the MP method is applied to actual 


chromatographic peaks. The assumptions made, and the mathematical 
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details of the peak model and the computation method used are presented 


in Appendix I. 
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EXPERIMENTAL 


Chromatograms were digitized by taking tnterval-area samples 
(defined in Chapter 1). An Aerograph Digital Integrator Model 471-42 
operated in the manual integrate mode was used for this purpose. The 
interval areas were printed on paper tape and manually transferred to 
punched cards. The remote readout terminals were connected to an 
Adjustable Time Delay (Potter and Brumfield, Model CHD) through a 


100-ohm resistor and a 0.02-MFD capacitor as in Figure 4. 


Adjustable 


Remote Time Delay 


Readout 
Terminals 
of Integrator 


Power Supply 


Figure 4. Schematic wiring diagram of circuit used to obtain interval- 


area samples of detector signal. 


The capacitor limited the integration dead time to less than 


60 usec as measured with an oscilloscope. A Gow Mac Model 9999-C 


Power Supply provided constant voltage for the Adjustable Time Delay. 


=] 
The sensitivity of the integrator was specified to be: 4.5 “10 


volt sec/count. The sampling interval was measured independently 


for each chromatographic run by dividing the elapsed time, as 
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indicated by the integrator's internal clock (line- frequency 
synchronous), by the total number of samples taken. Drift in the 


sampling interval over a 14-hr period was less than 0.1% 


The thermal conductivity cell was a Gow Mac Model 9285 
semi-diffusion pretzel cell with type-W tungsten filaments mounted 
ina TR IIB temperature-regulated case. A Gow Mac Model 9999-C 
Power Supply supplied bridge current. A two-diaphram pressure regula- 
tor supplied helium carrier gas at 50 psig. Sections of thermometer 
capillary along with a needle valve controlled the flow of helium 
through the sample and reference sides of the detector separately. 
The flows were adjusted to 40 ml/min by use of a soap-bubble 
flowmeter. Sample of 1 ul of liquid plus 1ul of air were injected 


on-column with a Hamilton 10-yp1l syringe having a Chaney Adapter. 


Most of the data were obtained with a 4-m, 0.3l-cm 
i.d., copper column packed with 10% SE-30 on 60 to 80 mesh 
Chromosorb P. The column temperature was held at 80.10°C to within 
an estimated 0.05°C in a Colora water-bath thermostat. The 
temperature of the thermal conductivity cell was set at approximately 
20°C with. a ae ae voltage transformer rather than the internal 


temperature controller. The pressure drop across the column was 


139 cm Hg, barometric pressure varied between 68.7 and 69.1 cm Hg, 


detector current was 250 ma, and the sampling interval was 


hehe! O24! 0.002 sec. 


The following digital filter (the necessity for filtering 


the data is discussed in Chapter 1) was used for both minimizing the 
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74 
sampling end-point errors and for whitening base-line noise for 


estimation of power-density spectra: 
h(2) = &(2) - [(6(2-1)+6(g+1)]/2 (88) 


where 6 represents a delta function and 2 represents reduced 

time in units of the sampling interval, At . The error the last 
data samples were given the value zero after filtering. The filtered 
data were corrected for both the filtering and the distortion caused 
by interval-area sampling by multiplying the Fourier transform of the 
filtered data by (ms)/[2 reyes wig , where § is reduced frequency. 
The corrected Fourier transform of the data at zero frequency was 
defined to be zero because the mean of the base-line noise, being 
arbitrary, causes the power-density of the noise at zero frequency 


to be extremely large. 
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POWER-DENSITY SPECTRA 


Curve A in Figure 5 is an estimate of the power-density 
spectrum of the base-line noise for the gas-chromatographic system 
described above. The estimate was obtained from five base-line 
sections of 512 interval-area samples each. The sections were 
taken between chromatographic runs over a period of 14 hr. The five 
Separate estimates of the power-density spectrum were averaged; the 
averaging process smooths the estimate and is equivalent to convolution 
with a Bartlett spectral window (34). A Hanning spectral window used 


for further smoothing is described by 
W(s) = 6(s)/2 + [6(s-1/N) +6 (s+1/N)]/4 (89) 


where W is the number of data samples. The power-density spectrum 
of the base-line noise after filtering (Equation 88), but before 
smoothing and correcting, was not flat, but decreased slightly at low 
frequencies. The filter did not give a completely flat spectrum, 

and smoothing this uncorrected spectrum caused a positive bias at low 
frequencies. Nevertheless, the error is not so large as the negative 
bias that would be introduced if the spectrum were smoothed after 
eorrecting tor filtering. ‘The ooey (total variance) of the noise was 
7 as and the qutocorrelation width 242 sec. The power of the 


2 
individual sections of base line ranged from 0.4 to SS ative gy atict 


the autocorrelation widths ranged from 115 to 243 sec. The variation 


in the power is a measure of the degree of nonstationarity in the 


system. The observed variation was much greater than would be 


expected of a stationary process. Nearly all the variation occurred 
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Figure 5. Power-density spectra of base-line noise from a thermal 


conductivity detector. The dashed lines indicate the Nyquist 
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ad 
pe the lowest. frequencies..(0,.001.. to. 0:004:Hz)<- Tests with computer- 
generated peaks (discussed in detail below) show that these large 
variations in power have little effect on MP estimates of peak 


parameters. 


Curve B in Figure 5 is the estimate of a power-density 
spectrum obtained for base-line noise with the automatic temperature 
controller of the detector in operation. The controller introduced 
an oscillation in the base line with an amplitude of about 10 uv 
peak to peak. The oscillation accounts for the slight shoulder on 
curve B at about 0.005 Hz, which corresponds to the cycle frequency 
of the controller. The power and the autocorrelation width of 
curve B is 110 ee and 280 sec. Curve B is less smooth than 
curve A in Figure 5 because it was derived from only one section of 


512 samples instead of five as for curve A. 


The base-line noise of the gas-chromatographic system 
(curve A in Figure 5) is characterized by intense noise at low 
frequencies. Above 0.02 Hz the noise is much less intense and is 
essentially white. The white part of the spectrum has an average 
power density of 0.017 ws /Hz . About 25% of the white noise 
component is caused by quantization which is a noninstrumental source 
of noise. The remaining instrumental white noise has an average 
power-density of 0.012 WOE ; the primary sources of this noise are 
probably the amplifier of the integrator and the power supply of oe 
detector. The relatively large proportion of quantization noise could 


be decreased by either amplifying the detector signal or reducing the 


quantization interval. 
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Interval-area sampling attenuates the noise at frequencies 
above the Nyquist sampling frequency [1/(2At) , where At is the 
sampling interval], but just above and less than twice the Nyquist 
frequency there is little attenuation (Figure 3, I). Sampling causes 
the power density of the noise in this region to overlap the power 
density just below the Nyquist frequency (Figure 3, VI). The overlap 
accounts for part of the upturn of spectrum A as the Nyquist frequency 
is approached, because the instrumental white-noise component probably 
extends through the region where overlap occurs. The remainder of the 
upturn is caused by the application of the correction for distortion 
due to interval-area sampling to quantization noise. Quantization 
noise is introduced after sampling and therefore does not require 
correction, but it cannot be separated from base-line noise, which 


does require correction. 


The random variation of the sampling interval (jitter) was 
measured by taking a number of sets of noise data with the zero- 
offset control of the detector bridge adjusted to give average outputs 
ranging from O to 30 mv . The power density, averaged over 
twenty points nearest the Nyquist frequency and nearest one half of 
the Nyquist frequency, when pletted against the square of the signal 
average, gave two straight lines with little scatter. The slopes of 
the lines were substituted in the equation for jitter noise (Equation 
41) and the equation for noise due to variation in integration dead 


time to give 40 psec as the standard deviation of jitter, and 


90 usec as the standard deviation of the variation in integration 


dead time. The result for jitter is reasonable and probably accurate, 
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but the result for variation in dead-time is larger than the integration 
dead-time itself. This violates the assumption made in the derivation 
of the effects of variation in integration dead time so the value 

90 uSec is not an accurate estimate of the true value. A result 
faas high could occur if the p.d.f. of the variations was non-normal, 
but more likely was caused by a noise source other than variation in 
integration dead time that also has a white power-density spectrum 
that increases with the amplitude of the detector signal. The noise 
from the integrator amplifier undoubtedly increases to some extent 
with increasing level of input signal. The power density of the noise 
at low frequencies also increased with signal level, but not so 
regularly and not at the same rate as at high frequencies. Further 
runs with a well insulated mercury cell and a voltage divider as a 
signal source showed conclusively that an important part of the power 
density of the noise at low frequencies was from the integrator 


amplifier (see Table VIII). 
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ESTIMATOR EVALUATION 


The estimation method was tested by adding computer- 
generated peaks with known parameter values to the sections of base- 
line noise used for estimating the power-density Spectra. Other than 
round-off errors in the computer-generated peaks (which were negligible) 
and numerical errors in the estimation process (see Appendix II), the 
base-line noise was the only source of error in the analysis of these 
simulated chromatograms. The computer-generated peaks were not 
interval-area sampled, hence the sections of base-line noise and the 
power-density spectra were not corrected for interval-area sampling. 
This means that the average characteristics of the base-line noise of 
the system being tested are slightly different from those of the 
experimental system. The difference is small (the correction for 
interval-area sampling is relatively small), however, and since the 
MP method is independent of the type of noise (as long as an accurate 
power-density spectrum is used), the results of the test will be 
valid for the experimental system. The retention times of the 
computer-generated peaks were taken from a set of pseudo-random real 
(nonintegral) numbers with a uniform p.d.f. The estimated parameter 
values were compared with the il RTE through the usevor 
Student's t . Bias was considered significant when the 90% 


confidence interval about the observed value of Student's t¢ did not 


include the value zero. The variances of the estimates were compared 


by use of the efficiency, which is the ratio of the predicted to 


observed variance. The predicted variance, obtained from the 


Cramér-Rao inequality (Equation 69), is the smallest possible variance 
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of any estimator. 


The predicted standard deviations in Figure 6 are square 
roots of the averages of the predicted variances for each parameter, 
calculated for each set of peaks. The relative standard deviation 
of area and the absolute standard deviations of the other parameters 
are approximately inversely proportional to peak area over the range 
oie to 10° counts. This range would correspond to about 0.0003 
to 0.3 ul of normal octane. For areas or amounts greater than this 
range the curves in Figure 6 should continue as straight lines. 


Therefore, when base-line noise is the only source of error, results 


as precise as desired could be obtained by increasing peak area. 


Table I gives the actual performance of the MP estimator 
for peaks similar to normal octane, Many of the efficiencies are 
Significantly greater than one. Theoretically the efficiency of an 
estimator cannot be greater than one, but an experimentally measured 
efficiency can be greater than one for two reasons. The first is 
the randomness of both the observed variance and the expected variance. 
The variance of the observed variance is approximately inversely 
proportional tothe number of experimental runs. Because the peak 
model is nonlinear, the predicted variance of a parameter depends 
on the values of all the parameters and will vary as the parameter 
estimates vary. When the total information about the parameters is 
reasonably large, however, the variation of the predicted variance 


should be negligible in comparison with the variation or variance of 


the observed variance. From the ten runs furnishing the data for peak 


area in Figure 6, the relative standard deviation of the predicted 
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Figure 6. Variation of predicted standard deviation with peak area. 


Except for area, the parameters of the peaks have values similar to 
those of normal octane in Table IV. Area is in units of integrator 
counts; nae units correspond to about 0.03%. -of-a~-l-ul-—sample. 
The standard deviations of area, ee , are relative, and those of 

t= “and ‘0 are in Units of the sampling interval. The power— 


0) 


density spectrum is that of curve A in Figure 5. 
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TABLE I 


Effect of Peak Area on MP Estimator Performance® 


Briiciency (¢) 


Area 102 10 10> wee 

A 4.4 1.6 a 
eo 567 AGP 3.0 Aad 

ie} oie 2 IG 943 2.0 
K, eu! eyes iN) 204 
K, shee 5.0 2.9 2.5 
ss 6.2 2.9 1s DAES | 
Se G3 TON 5.8 238 


Cc 
Bias (Student's ¢) 


5 6 
Area 10° 10° 10 10 
A = a5 201 Le? TO 
£ 120 Oy OTL 0.6 
0 
Oo Sree 0.5 1.3 14 
K, -0.6 =O Or? 0.8 
K -O.1 1220 14 Jog | 
4 
Ke -0.8 onl Ge 0.5 
Ke AS 1.6 2.0 2.6 


“ren computer generated peaks were evaluated for each area. 


The peak parameters are as in Figure 6. 


rhe 90% confidence interval for ec is 0.53€, < el < 2.7€ os 


where Ey is the observed efficiency. 


“the 90% confidence interval is ¢+1.8. 
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variance was 10% at igs and 0.01%. at 10° area units. The 
approximate relative standard deviation for an observed variance, on 
the other hand, is /2/n where ” is the number or runs, or 45% 
for ten runs. Confidence intervals for the estimated efiveiency, 
which is a ratio of variances, were calculated by use of Fisher's 

F distribution. The predicted variance in the numerator of the 
efficiency ratio was assumed to have an infinite number of degrees of 


freedom. 


The second cause of an estimated efficiency greater than 
one is inaccuracy in the power-density spectrum. For example, if 
the spectrum is mutliplied by a constant, the predicted variance, 
and therefore the efficiency, will be proportional to the constant; 
the observed variance, however, will not change. The positive 
bias at low frequencies of the estimate of the power-density spectrum 
Probably ast theimajor cause of the efficiencies in Table I being 
greater than one. The error in the eee ee spectrum should not 
have serious consequences; it should have little effect on bias of 
the parameters (except possibly for extremely broad peaks, discussed 
below) and merely make the predicted standard deviations slightly 


conservative. The efficiencies still can be used to detect trends 


in the performance of the MP estimator. 


Effect of Peak Area 


Contrary to’ the theoretical prediction in Chapter 1, Table I 


shows that as peak area and information increased, efficiency seemed 


to decrease and bias to increase. (The amount of information about a 
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parameter is roughly equal to the inverse Square of the predicted 
standard deviation.) In view of the randomness of efficiency and bias 
these variations may be coincidental. The steady change in bias from 
negative to positive, and the less pronounced decrease in efficiency 
with increasing area is probably the result of numerical errors in the 
estimation process rather than a property of the MP estimator itself. 
As the amount of information increases, the demand for numerical 
eecuracy also increases. But all the data in Tables I, II, and III 
were obtained with approximately the same numerical accuracy in that 
the probability of the estimates was close to the probability of the 
most probable values of the parameters (same termination constant, 

see Appendix II). A better procedure might be to vary the termination 
constant with the amount of information so the numerical accuracy 


can increase as the information increases. 


Effect of Peak Width 

Again contrary to the predicted trend, Table II shows that 
efficiency generally seems to decrease as peak width decreases and 
as information increases. The cause of the disagreement with theory 


is probably the positive bias in the estimate of the power-density 
Spectrum at low frequencies. As peak width increases, the information 


for all the parameters concentrates at low frequencies. Thus, the 
information matrix, and consequently the predicted variances, for wide 


peaks will be influenced more than for narrow peaks by the bias in the 


power-density spectrum. The decrease in apparent efficiency with 


increasing peak width is caused by a determinate error in the predicted 
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efficiency is caused by a determinate error in the observed variances. 


As the total amount of information about a peak decreases, 
the posterior p.d.f. flattens and eventually more than one maximum 
May appear. The base-line noise may form peak-like curves, which will 
appear aS maxima in the posterior p.d.f. and will be especially 
evident along the retention-time axis. Predicted variances, which 
depend on the curvature of the posterior p.d.f. near the absolute 
maximum, are accurate only when the posterior p.d.f. is normal. 

A posterior p.d.f. with multiple maxima is definitely not normal, 
and the actual variances will be greater than predicted. Therefore, 
as information decreases and the detection limit is approached, 
efficiency should decrease, though the decrease does not begin within 


the ranges of area and peak width in Tables I and II. 


Although the bias of most of the parameters in Table I is 
insignificant, there seems to be a general trend of increasing bias 
with increasing order of the parameters, which might be expected 
because of the increasing nonlinearity of the high-order parameters. 
If the test for bias were based on the predicted rather than the 
observed variance, the bias of all the parameters would be significant 
In Table II the 


because the efficiencies are all greater than one. 


bias appears to be independent of the peak width. 


Number of Data Points 
The effect of the number of data points used in the 


estimation on bias and efficiency is shown in Table III. Evidently, 


the rule of thumb derived in Chapter 1 (the number of data points that 
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AMANB EEO! RICE 
Effect of Number of Data Points on MP Estimator” 
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Grive runs were made for each value of N (number of data 


points used in the estimation). Same conditions as for Table II. 


2 is the autocorrelation width of 


= here I 
T’. 80 + T. , whe a 


the noise. 


&8 


Od hE tat 


jm 
C 
* 
Lo) 


&%6D 2to xsdmpn) Vi 


Qo «bite noite lox 1opOIUs enh thee Rs 


et ist: 


me | £.0- 2.9 


Jo anivy” asks, cae 48 
tO? es snot tkbnp3, eae 


| eh See ; 


tt At B24 ste. 28: saa Re 


ST = ee ipa = & 


= 


ee — ee erent + — : ae [= A ot Ti 
I . a 
i) 
_ ae ese 5 
oe a 
ty el) 
} Ef Ko. 
i HH; oO 
ry 7 APU 
-: 
s. 
} = ae 
Y+U en. Ah 
0 ras © 
> €,.0 e.G 
$ mb ,£00..G 
f£5 2 


620.0 
j ei 


89 
cover eight peak standard deviations plus the autocorrelation width of 
the noise)gives an adequate number of points for the noise studied. 
Half that number decreases the efficiency for any of the parameters by 
no more than 50%. Bias seems to stay constant as the number of 
points decreases, but since the efficiency decreases (because the 
observed variance increases), the value of the bias in units of the 
parameters actually increases. The swum of wetghted squared devtations, 
represented by the symbol @ in Table III, should theoretically be 
nine less than the number of points, W (see Appendix I). A constant 
difference between @ and WN as W decreases could be explained by 
bias in the estimate of the power-density spectrum. The actual value 
of @ decreases more slowly than W in Table III, probably because 
of end effects arising from the chromatographic signal being nonzero 


near the ends of the interval over which samples were taken. 


Effect of Nonstationary Noise 

Nonstationarity was examined by analyzing sets of five to 
ten computer-generated peaks added to the lowest- and highest-power 
base-line sections used to obtain curve A in Figure II. The two 
sections were taken about six hours apart. The efficiencies were 
usually greater by a factor of two on the low-power noise than on the 
high-power noise for peak widths ranging from ten to twenty sampling- 
interval units, but there was no consistent difference in bias. 
100 


Even though the power of the high-noise section was more than 


times the power of the low-noise section, the values of @ were not 


greatly different: 520 and 420. The power of noise, therefore, 


is a poor indicator of the effect of the noise on parameter estimates. 
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Moreover, the nonstationarity encountered did not seriously affect 


the estimation process. 
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CHROMATOGRAPHIC PEAKS 


In this section results obtained from actual chromatographic 
peaks are presented. No effort was made to optimize the performance 
of the chromatographic system. The column was designed only to give 
peaks broad enough to minimize sampling errors with the integrator- 
sampling system operated at the maximum practical sampling rate. 
Tables IV, V, and VI _ summarize the results of measurements of two 
hydrocarbon mixtures. One mixture contained about 3.6% by weight 
of each of normal heptance, 2,3,4-trimethylpentane, and normal 
octane. The other mixture contained about 4% of each of normal 


heptane, toluene, and normal octane. The solvent was 2-methylpentane. 


TABLE IV 


Means of Peak Parameters 


a lene Cc K 

= eg " K Ky os 6 
ae 1.025 262.6. 3.89 0.107 0408 @ 0Ml5 @ 0.74 
heptane 
pecmat hot 6.54579 O 7.24 0.060 “80404 m@ los 2008 
octane 
“0: aaa 919 391.7 5.42 0.045 0404. Omlo = 0,08 
methylpentane 
Toluene .919 A13—5 5.88 0.464 0.40 Oe2 1 0. 52 


ee 


“Relative area response per gram. . 
Areas were normalized, divided by the weight of the component, 


then averaged (not relative to any particular compound). 


Det retention time. 


“th units of the sampling interval. 
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Standard Deviations 


The observed standard deviations in Table V are all greater 
than those predicted. The low efficiency could not have been caused 
by the efficiency inherent in the MP estimator because the efficiency 
obtained for computer-generated peaks was generally greater than one. 
The predicted standard deviations account for the effect of base-line 
noise only; the low efficiency must be the result of non-base-line 
noise caused by fluctuations in flow rate and pressure of the carrier 
Gas, temperature of the detector cell and the column, etc. 

The observed and predicted standard deviations for relative 
weamiarhzed area agreed with a factor of 3 to 6. In spite of 
the crude chromatographic system, the standard deviations were low by 
usual criteria. With an apparatus designed to control non-base-line 
noise, Goedert and Guiochon (12) attained a relative standard 
@evratron of only 0.6 x tom for a peak with roughly the same width 
but twenty times the area, in absolute units, of the normal heptane 
peak described in Table IV. The predicted standard deviations for 
relative normalized area were calculated from the predicted standard 
deviations of the individual area measurements”. 


*het Ay pees area of the component of interest and 


AyA,,---+1A be the areas of the other components. The normalized 
n 
area is 
= (90) 
Ay As 
ie . . 
where S$ = ) A. . The square of the standard deviation of the 
ga 
normalized area is 
lay ote 
(nat)? = [(Staa,) “4a YM) 1/8 (91) 
dL j=? 
a Je . . 
Where S§' = ) yy rand: mild eee standard deviation of component 
. J J 
ee 


tee The areas are assumed to-be statistically independent and this will 


apply when the peaks are far enough apart so the information functions 
for areas do not overlap in the time domain. 
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The observed relative standard deviations are greater for absolute 
-3 
area (34 < 10 for normal heptane) than for normalized area because 


of the sample-injection technique. 


The observed standard deviations for retention time in 
Table V are about 1000 times those predicted, probably because of 
low-frequency fluctuations of the carrier-gas flow rate and the column 
temperature. Low-frequency fluctuations should not affect relative 
retention-time or retention-index measurements as much as absolute 
retention-time measurements. The mean, observed standard deviation, 
and predicted standard deviation for the retention index of 2,3,4- 
trimethylpentane were 754.61, 0.022, and 0.0019; for toluene these 
were 762.66, 0.026, and 0.0018. The observed standard deviations 
for retention index were indeed only ten times the predicted values. 
Precision estimates are seldom included in published retention-index 
data. Retention indices are usually reported to the nearest whole 
number and seldom to the nearest tenth (34). Precise measurements of 
retention indices may not be necessary because absolute retention 
indices are difficult to reproduce on different columns. Nevertheless, 
precise values are needed to compare the retention indices of different 
compounds on the same column and to measure the rate of change of 
retention index with column temperature. The highly precise estimates 


of retention time obtainable by the method presented here should be 


valuable for these purposes. 


The predicted standard deviations for the retention indices 


reported above include the contribution from the predicted standard 
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deviation of the retention time of the air peak*. Although the air 
peaks were too narrow (o equal to 1.0 sampling-interval units) to 
obtain a good fit, the MP estimates of the air-peak retention times 
were used in the retention-index calculations. The derivative with 
respect to retention time of the function M, which is used to obtain 
the MP estimates (see Appendix II), is antisymmetric in the 


frequency domain, whereas the overlap error in the frequency domain 


* 
BOteme. sot Re ere be the retention times, and At, 
a n Be n+l a 


Avs At. At ed be the standard deviations of the retention times of 


@impoche «7*th normal ‘hydrocarbon, compound «x ., and thee @+l'th 
normal hydrocarbon. The retention index of compound x is then 


NE = 100” + 100N/D (92) 
ARB 
where 
N= log(t_~- t_) - LOGIE ots) 
and 
- - -¢t 
Die log(t, 5 t.) log (t, =) 


The square of the standard deviation of the retention index is 


2 100° are u oS a 2 
(Gee eee ets eas le ee es a ace 
eS x a N+1 “a 
(N-D) At Ss 2 
prea of ie, [pe Dish. coeetion iets (93) 
# -t nT. e -¢t t -t a 
hen x a n+l “a na 


Equation 93 gives the minimum variance derived from the Cramér- 
Rao inequality (Equation 69). Equation 92 is the Ges neon or a 
retention index, but the direct substitution of RASA SAREE) eh estimates 
into the defining equation does not give an MP SE Beh of the 
retention index. Similarly, the definition of peak area is unrelated 
to the MP estimator of area. These remarks also apply Foo the 
estimates of relative normalized area discussed in the (CLG ASS footnote, 
that is, the estimator of normalized area based on the definition of 


normalized area may not be the best estimator. 
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caused by the peak being too narrow is symmetric and should have little 
effect on the retention-time estimates of the air peaks. The error 
curves (the difference between observed and fitted peak, see 
Appendix I) for air peaks were indeed symmetric about the retention 
times, showing that the error in the fit was primarily in the area and 
width parameters. In any case, the effect of small errors on the 
retention time of the air peak is negligible in the retention index of 


peaks far from the air peak. 


Base-line noise accounted for about 10% of the observed 
standard deviations for the peak width and the cumulant measurements. 
Part of the remaining variation may be due to noise that increased 
with signal amplitude; this noise would cause variations in peak 
Shape that would not affect area and retention time so much as it 


would the higher-order parameters. 


In Table VI, the observed correlations between parameter 
deviations from the mean were calculated by use of the Olkin and 
Pratt correction for bias (35). Most of the observed correlations 
agree with the predicted correlations within the 90% confidence 
interval. The correlations between errors in retention time and peak 


width, however, are consistently and significantly higher than the 


predicted values of zero. The high correlations are caused by the 


strong physical relation between retention time and peak width that 


was not included in the peak model. If this relation were known 


exactly, the retention-time and peak-width parameters could be 


replaced with a single parameter. 
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Low-frequency fluctuations in the flow rate of the carrier 

gas with a period longer than peak retention times would introduce 
a proportional error into the retention times of all the peaks in 
a chromatogram. The high positive correlation (greater than 0.97) 
found between deviations of retention times for all pairs of peaks on 
the same chromatogram confirmed the effect of extremely low-frequency 
flow fluctuations. Peak-width errors were affected less strongly; 
the correlations between errors in the widths of peaks on the same 
chromatogram were greater than 0.56 for all pairs of peaks. 
Cumulants were unaffected; the correlations between errors in the 


cumulants of different peaks were near zero. 


The sum of weighted squared deviations, @ , for the peaks 
described in Table IV exceeded the expected values by 10 to 1000 
times. Part of the excess deviation is due to the nonstationarity of 
Pifemnoise; the values of @ varied by a factor of 10 from run to 
run for the same component. The two major contributions to the large 
values of @ are incomplete agreement of the true peak shape with the 
peak model (poor fit) and interfering peaks (impurities that give 
small peaks near the peak of the principal compound). Interfering 


peaks will be discussed first. 


Interfering Peaks 
A typical error curve for toluene (Figure 7, I) shows what 


appears to be a single interfering peak with.a maximum...3./0,7.greater 


than that of the principal peak. The effects of interfering peaks on 


MP estimation were studied by adding principal and minor pairs of 


computer-generated peaks with similar parameters (except for area and 
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Figure 7. Error curves of chromatographic peaks. I, toluene; II, 
2,3,4-trimethylpentane; III and IV, normal octane; V, normal 
heptane; VI, normal nonane with severe drift. The number at the 

top of each figure is the approximate height of the peak in microvolts. 
The dashed lines are the profiles of the peaks. The mean values of 

the parameters of the peaks in I, II, III, and V are given in 


Table IV. 
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retention time) to a section of base-line noise. Error curves 
Similar to Figure 7, I were obtained with peaks comparable in size 
and shape to the toluene pair. The height of the peak on the error 
curve was found to be proportional to the true area of the interfering 
peak. Relative peak heights give the area of the interfering toluene 
peak, if it is a peak, as 2.7% of the area of the toluene peak. 
Similarly, the error curve for 2,3,4-trimethylpentane (Figure 7, II) 
shows an interfering peak that would be about 0.5% of the principal 
peak. The normal heptane peak appears to have two interfering 


Gomponents (Figure 7, V) of about 0.3%: .and 0.4% . 


The effect of an interfering peak on @ can be calculated 


from the amount of information in the data about area (I defined 


rei SY 
by Equation A-7 in the Appendix). The amount of information for area 
is the sum of weighted squared deviations of a peak with unit area 
from a constant zero base line. Therefore, the amount of area 
information for the toluene peak times the square of the area of the 
interfering peak should give a crude estimate of the contribution of 
@ of the interfering peak. The median value of @ for toluene was 
4.4 x 104 , and the estimated contribution of the interfering peak is 
450 x 10° . The values of @Q obtained for computer-generated 
interfering peaks agreed with the estimated value within an order of 
magnitude. The agreement was not so good for heptane and 2,3,4- 


trimethylpentane, however. For these peaks the interfering peaks 


accounted for only 10% of the excess in the values of @. 


The detectability of interfering peaks was studied by running 


mixtures of 2,2,4-trimethylpentane (principal) and heptane (minor) 
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diluted with 2,3,4-trimethylpentane. The concentration of 2,2,4- 
trimethylpentane was 3.6% by volume; the separation between the 
principal and minor components was (3.0 + 0.4)qo (resolution 

0.7 + 0.1). The left part of Figure 8 shows the peaks as they 
appeared on a strip chart recorder. Where the concentration of 
normal heptane is 10% of that of 2,2,4-trimethylpentane, heptane 
appears as a Slight elongation of the tail of the principal peak. Below 
10% there is no visual indication of an interfering peak. With the 
use of error curves, however, an interfering peak is evident at 0.13 
(not shown in Figure 8), and at 0.6% there is no doubt that an 
interfering peak is present. The use of error curves increases the 
detectability of an interfering peak by a factor of at least 10 in 


this case, 


Relying on the detection of interfering peaks in the manner 
described above may be hazardous. If the peak model is too general, 
that is, /if the parameters contain insufficient prior information, 
the peak model may fit multiple peaks as a single entity. A study 
of computer-generated peaks using up to the sixth cumulant in the peak 
model, revealed that an interfering peak could be fitted along with 
the principal peak (incorporated within the model) without leaving a 
trace of the interfering peak on the error curve. The factors 
Me senting the incorporation of an interfering peak are the relative 
areas of the peaks, the difference between their retention times, 
and the total amount of information available about the principal 
peak. The information in the data about the principal peak depends on 


the shape and size of the peak relative to the power-density spectrum 
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10% 


2.8% 


Figure 8. Detection of an interfering peak by use of chart recordings 
(left) and error curves (right). The percentage figures refer to the 


concentration of the interfering peak relative to the principal peak. 
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of the base-line noise. For example, the information about a peak 
decreases with increasing peak width. The interfering peaks ina 

set of computer-generated peaks were visible on the error curves when 
the peak widths were equal to 2 sampling-interval units, but were 
incorporated within the peak model and did not appear on the error 
curves when the peak widths were equal to 10 units. (The area of 
the principal peak was 10° integrator counts, the Separation was 

30 (resolution 0.75), and the area of the minor peak was 0.13% 

of the area of the principal peak.) Although an error curve may not 
indicate an interfering peak when it is incorporated within the peak 
model, a peak that is incorporated will have a distinct effect on the 
values obtained for the cumulants. If the cumulants of the pure 
principal compound were known, an interfering peak could be detected 
by comparing the cumulants of the pure and the impure compounds. It 
would be just as well, though, to introduce the prior information 
about the pure compound directly into the peak model or into the 


Piror 0.d,f,., of the cumulants. 


When error curves of computer-generated peaks have the 
appearance of the error curve of toluene in Figure 7, I, all the 
parameter estimates have positive bias except area, which is slightly 
less than the true value. When the retention time of the interfering 


peak is less than that of the principal peak, giving an error curve 
the mirror image of Figure 7, I , the effect of the interfering peak 


on the even-order parameters is the same as when the retention time 


of the interfering peak is greater than that of the principal peak, 


whereas the odd-order parameters have negative bias. The values of 


the cumulants obtained for toluene are significantly greater than those 
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of 2,3,4-trimethylpentane (Table IV). If the peak in the error 
curve for toluene were due to an interfering peak, the difference 
between the cumulants of toluene and 2,3,4-trimethylpentane could be 
attributed to the bias in the estimates of the cumulants of toluene 


due to the interfering peak. 


Limitations of Peak Model 

The above results on the detectability of interfering peaks 
are self-consistent and valid with respect to the a prtort information. 
Nevertheless, attempts to isolate an impurity in toluene by fractional 
trapping were unsuccessful. The toluene did not contain an impurity 
with a concentration as high as 2.7% . Furthermore, error curves 
obtained for benzene were similar to those for toluene, and the benzene 
would be unlikely to have an impurity with the same relative area and 
Peeenition as ah impurity in toluene. The peaks on the error curves of 
benzene, toluene, and possibly normal heptane and 2,3,4-trimethylpen- 
tane were due to residuals that the peak model could not fit. These 
peaks could be fitted if the peak model were extended to include more 
adjustable parameters, but a better approach would be to use a different 
basis function. A Poisson basis function might be suitable for the 


peaks of the aromatic compounds studied here. 


The cumulants of a Poisson function are all equal toa 
positive constant. Therefore, attempting to fit a Poisson peak with 
a Gaussian basis function and a finite number of cumulants as adjustable 


parameters will introduce negative bias into the nonadjustable 


cumulants, which are assumed to be Zero. This negative bias will 


Propagate downward to the low-order adjustable cumulants through the 
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correlation between errors in high- and low-order Pe (An 
interfering peak, in comparison, will cause a positive bias in ne 
low-order adjustable cumulants.) Therefore, the true values of the 
cumulants for toluene may be Significantly greater than the values 

in Table IV. Even though they are biased, the difference between the 
estimated values of the cumulants of 2,3,4-trimethylpentane and 
toluene are clearly evident. The numerical differences between these 
peak-shape parameters show that cumulants can be used to aid 
identification. The visual difference between the shapes of 2,3,4- 
trimethylpentane and toluene on a strip-chart recorder was relatively 


small. 


The error curves obtained for normal octane (Figure 7, III 
is a typical example) indicate a poor fit without a straight-forward 
explanation. Estimating the values of the cumulants beyond the 
sixth order should improve the fit. Without actually using an 
expanded peak model, the effect of including additional parameters 
was studied by examining the error curves of computer-generated peaks 
with nonzero values of cumulants up to the tenth order. The error 
Curve in Figure 7, III is characterized by an envelope of high- 


frequency oscillations, the major component of which has a distance 


between nodal points of 0.50 . The error curves obtained when the 


peak model was terminated at the next lower cumulant were similar to 


: ‘ he t inal 
the curves in Figure 9 for K, to Ke As the order of the termin 


cumulant increases, the distance between nodal points on the error 


curve near the retention time decreases, but only from Lede = 2% K, 


eo L005 for. Hence, if the error curve of normal octane was 
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Figure 9. Additive effect of cumulants on a Gaussian peak in the 


time domain. Each curve is the difference between a standardized 


peak shape function in the form of Edgeworth's series with one non- 


zero cumulant and a standardized Gaussian peak shape function. The 


values of the standardized cumulants, K, to Ke praLe: sO, 
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caused by use of an incomplete model, the normal octane peak would 
require a model with adjustable cumulants up to perhaps the twentieth 
order. Such high-order cumulants could not reasonably cause 
oscillations of the magnitude seen in Figure 7, III because of the 
decreasing amount of information in a cumulant, and hence the decreasing 
influence of a cumulant, as order increases. A number of cumulants 
considered as a group may be significant, however, as in the case of 

a Poisson peak where each cumulant is equal to the same constant. If 
there was a functional relation between the cumulants of a peak, the 
parameters of the function could be estimated rather than a large 


number of individual cumulants. 


As peak area and information decreases, the influence of the 
nonadjustable cumulants decreases, and better fits can be obtained. 
The error curve for octane (Figure 7, IV) at one tenth the area of the 
peak Giving the error curve in Figure 7, III indicates a perfect fit. 
At lower areas, however, the estimates of the cumulants are much less 


precise and are valueless for the purpose of comparing peak shapes. 


Error curves similar to Figure 7, III were obtained when 
white noise with power proportional to the square of the signal 
amplitude was added to computer-generated peaks in a simulation of 
jitter noise. But to produce oscillations of the same magnitude as 


those in Figure 7, III the standard deviation of jitter would have to 


be more than ten times the value found from the study of power-density 


spectra. Furthermore, the major features of error curves of normal 


octane were reproducible from run to run; reproducibility would not be 


possible if the oscillations were random. 
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CONCLUSION 


In general, if enough data points are used, if the base-line 
noise is reasonably stationary, and if all the nonadjustable parameters 
are either zero or insignificant, the MP estimates of the first 
seven parameters of Edgeworth's series are unbiased, and the variances 
of the parameters due to base-line noise are approximately the same as 
the predicted variances. The method is insensitive to severe low- 
Pecguency noise or drift; Figure 7, VI shows a perfect fit in the 
presence of an intense low-frequency oscillation caused by the 


temperature controller of the detector. 


The most important limitation of the method is that the 
conditions listed in the above paragraph must be known to apply to the 
system under study before the estimation is carried out. All 
measurement methods are subject to similar a prtor?t conditions, but 
the MP method brings these assumptions to the foreground where they 


can be clearly defined and tested. 


One of the most worrisome assumptions is that of stationarity 
of the base-line noise. We found the greatest effect of nonstationarity 
to be in the accuracy of predicted standard deviations. If more 
accurate standard deviations were required, the effect of nonstation- 
eraityeand the:effect; of errorsin) the measurement of power-density 
spectra could be reduced by approximating the power-density spectrum 
with a simple function of one or two parameters. These parameters 


could then be estimated along with the peak parameters to give an 


accurate estimate of the current power-density spectrum. 


7 
ioave © 


edt ‘eae Vor, | 


® 
yr , 
iii-sfed ad 3h \bsew Ste =301204 sine dpdans ai 
’ - 
exstoncved olds teppteaan, ade far 2). be 
“— 18) pen patch OM 
= acts iin ots 
= BU LRGs “sS 
Ste A 


tamsand et eit sti Aut exonatiane 
rag a ekona’ Iv: oe Hsiaie aie ae 
ioteel bisso isaeupext-meh obebbet. 03 


a & =a 


> 
rane o ‘ona $6 ‘ideretee” xt 
ot ia 


~ 


~*~, Di wh 


wit ects: div forte ett te ieee te FLae Sf — feom ¢ 
“ 
~~ al 


: i Am 
of3 of yYloqus of “nworntt Se stabs ge? Asay ever att Rant 
pages 


4 peas > 
» 


‘ri . 


a be 
‘cd DALY AD St aon tite ‘gm ‘eit a] oa 
a ae 


‘We ry " : ’ Sern ‘a can nd . 
Sud, enor Lbnca: Beales vst came om gne tide it horton 
< " 


AS 
: VSP ee 
mt 


yer? sus briiotpstoipeds as ceoesuiuel Sep, 


Weissnoiseys Jo seid Be znd ene 3 ie seh saat +a 
| . Ci 
veitenoiszsetecc: in Soskts sess ists feive a pode 


eon 4I .snokseteeb ohustivtecte hess eben te ‘ 
Le 


fa 


“weddsienon I. = sits Siobaiies otew sven bees 
yiderob “toyed 3 2 Acuna ‘ery ot sai % deste as 
arts 2s yF Lecuiti-ewag ae arta ee bean 9 | £ 
cesta S2out aan ons , Se ea: > nie yi a 


Gi 
1d aed Beas od 


Kereyete 


aha ie 


, 


Te 


110 


There is some reluctance in the literature about measuring 
high-order moments (8,9) because of the increasing influence of 
base-line noise on measurements of moments of increasing order. 

Using the MP method to measure cumulants (a cumulant is equal to 
the moment of the same order plus a combination of moments of lower 
orders), however, we have shown that the decrease in precision is not 
so rapid as to make the measurement of the sixth-order cumulant 
valueless. Depending on the power-density spectrum of the base-line 
noise, the precision may in fact increase with increasing order (in 


Table Vs. for Ky and » Kan). 


. 


The MP method is extremely sensitive to the presence of 
small peaks near the principal peak. This sensitivity is a disadvan- 
tage in that measurements of cumulants for their own sake requires 
highly pure samples. Nevertheless, the ability to detect a seemingly 
invisible peak on the shoulder of a large peak is valuable in itself. 
It would be better, both for the purpose of measuring the shape of a 
peak and of detecting interfering peaks, to extend the peak model to 
include multiple peaks each with a separate set of parameters. The 
fitting of multiple peaks has been done using a Gaussian peak model 
(36) and a model based on Sa ialiesners a previously run peak (15) . 
Greater effective resolution and lower detection limits should be 
attained, however, if a more complete and accurately defined peak 


model were used, and if the characteristics of the base-line noise were 


taken into consideration. 


The data presented is for peaks for which there was no prior 


information about any of the parameters. Prior information about any 
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one parameter will decrease the standard deviation for that parameter 
and any other parameter that is correlated with it. For example, con- 
Sider the effect of complete prior information (zero prior standard 
deviation) about high-order parameters on the predicted standard 
deviation for the area of a normal octane peak. The relative standard 
deviations (x1000) are 0.82 (with KL and higher well known), 
o-6lr uwith Ke and higher well known), 0.61 (with Ky s" Sten), 
0.43 (A), 0.42 (A), Om234 (0), 0.23 (Co)s and 0.0 (4, and 
all other parameters well known). The most significant parameters 


for area are K_ ,K 


6 and -O .© This is) also indicated iby the high 


A! 
correlations between these parameters and area shown in Table VI. 

Prior knowledge of the shape and width of the normal octane peak would 
reduce the standard deviation for area by a factor of 3.5 . On the 
other hand, prior knowledge of the retention time would have no 

effect on the precision of area estimates. It must be emphasized, 
however, that these results are the standard deviations due to base- 
line noise alone predicted for an efficient estimator. Although the 

MP estimator is nearly efficient when the amount of information is 
large, it is inefficient near the detection limit where the information 


is small. Near the detection limit, prior information about retention 


time may be important for the estimation of area. 
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APPENDIX I 


MODELS 


The posterior p.d.f. of the set of parameters, a 4 given 


the data, zo is given by Bayes' theorem 


(a|¢) = ————_____ (A-1) 


where P(e) PSCne PrlOr. p.d.fs, OL the ,paramerers, |e) 
v 


Dee 
zig 
the p.d.f. of the data given the parameters, and P,. (&) a 
~ 
scaling or normalizing factor. Since the maxima of a function and the 
maxima of the logarithm of the function occur together, finding the 


maximum of a posterior p.d.f. is equivalent to finding the minimum 


of the function 


M(a) = srlicppst sion + constant (A=9) 


where the constant is independent of the true values of the parameters, 


x 

A priori information about the experimental system is 
embodied in the functional form of the p.d.f.'s in Equation A-l. In 
the chromatographic system used here, the a prior? information consists 
of the following statements. 


(a) Since there was no prior information about the parameters, the 


prior p.d.f. of the parameters, p (@) , is assumed to be uniform for 


each parameter. 


(b) The detector signal is the sum of the base-line noise and the peak 
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function, and any noise with an amplitude that increases with increasing 
signal amplitude is either absent or negligible. 

(¢). The base-line noise is a stationary random process with a 
multivariate normal p.d.f. and a known power-density spectrum, CY) : 
(d) The peaks being studied are accurately represented by Edgeworth's 
form of the Gram-Charlier type-A series, and all cumulants above the 
sixth are equal to zero. An accurate representation means that the 
difference between the true peak function and the estimated peak 


function is indistinguishable from the base-line noise. 


With this prior information the M function becomes 


baeipraoaals 


cv) dv (A-3) 


where X(v) is the Fourier transform of the chromatographic data and 


F(v) the Fourier transform of the peak function: 


2 
F(v) = Al(cos@ + t sine) exp[-(2r0v) /2 


(A-4) 
4 6 
+ K, (2t0v) i242 K, (2r0v) /720) 
where 
3 5 
cic - (2nt v) - K , (2rov) joes K, (2nov) /120 
z is retention time, o is 


Vv is frequency (Hz), A is area, 


width, and K3 to kK. are standardized cumulants. The elements of 


the information matrix are given by Equation A-5. 
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where t indicates the complex conjugate. Predicted standard 
deviations are defined as the square root of the diagonal elements 

of the inverse of the information matrix. The M function can be 
evaluated by use of a digital computer provided that, in the digitiza- 
tion of the signal from the detector, the sampling interval is small 
enough, and that the sampling time is long enough that the integral 

in Equation A-3 can be replaced by a summation over a finite 


frequency range. 


Terminating Edgeworth's series at the sixth cumulant in 
Equation A-4 limits the range of peak shapes that this model can 
eccurately represent. One of the conditions on a peak model, if it 
is to be a realistic representation of an actual peak, is that the 
absolute value of the Fourier transform of the function at all 
frequencies must be less than or equal to the value of the function at 
zero frequency (37) . In the peak model defined by Equation A-4, 
this condition ieee the value of KKK must be less than 1.6. 
When this condition is met, a plot of the absolute value of F(v) 
against frequency has a peak-like profile with a maximum at zero 
Frequency. If KK is greater than 1.2 , however, the peak will 
have secondary maxima. In the critical region where KK is less 
than 1.6 but greater than 1.2 , the shape of the absolute value of 
F(v) changes rapidly with small changes in either K, or K.. : 
Figure 9 shows the time-domain difference between peaks with the 


indicated values of individual cumulants and a Gaussian peak with 
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cumulants that are all equal to zero. As this figure indicates, the 
effect of the even cumulants on the shape of a peak is symmetric 
about the mean, and the effect of the odd cumulants is not exactly 


but nearly antisymmetric. 


From Equation A-3 , when the qa prtort assumptions are true, 
the minimum of M is the sum of a number of random variables, each 
of which is normal with unit variance. Twice the minimum value of M, 
represented by the symbol @ , has a chi-square distribution with 
degrees of freedom equal to the number of data points less the sum 
of the number of parameters being estimated (seven) and the number of 
degrees of freedom lost by filtering (two, when the filter described 
by Equation 88 is used). The mean or expected value of @ is equal 
to the number of degrees of freedom, and if Q is significantly 
greater or less than expected, the a prtort assumptions are unlikely 
Rompe true: A graph of the error curve [x (t)-f(t,@)1, where @ 
represents the estimated values of the set of parameters, can be used 
as a qualitative test of the peak-model assumption. If the peak model 
is correct, the error curve in the region of a peak should be 


indistinguishable from the base-line noise before and after the peak. 


The curves in Figure 9 are error curves that would be obtained if 


there were no base-line noise and if individual cumulant estimates were 


in error. 
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APPENDIX IIL 


MINIMIZATION 


In this study the process of finding the minimum of the 
M-function in Equation A-3 was started by finding crude initial 
estimates of the parameters. The initial estimates of the cumulants 
were zero. The approximate area, retention time, and width of the 
peak were determined by use of the definitions in Chapter 1 after 
subtracting a linear base line from under the peak. The end points 
of this base line were taken as the center of a five-point moving 
average when, going downward on each side of the peak, the difference 
between a data point and the preceeding moving average was less than 
a preset value. After the initial estimates were obtained, the set 
of data samples containing the peak was filtered, Fourier transformed, 
and corrected for the distortion caused by filtering and interval-area 
sampling. Then a parameter correction vector, A , was calculated as 


outlined below, and the previous estimates were replaced by new ones: 


= (A-6) 
ee ee 


This procedure was repeated until a minimum of M was reached. 


As a function of two parameters at a time, M makes a 
surface that can be pictured as a distorted bowl. A Como of 
constant height on this bowl would look like an elongated and distorted 
ellipse. The elongation was reduced and the accuracy of the computer 
calculations was improved by transforming the parameters into a new 


coordinate system in which the amount of information about each 
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parameter is unity (38) : 


, : = ° ag o 8 “lf2 oa 
oe a+ ( ip (A-7) 
Be -1/2 : 


A ball released on the M-surface would begin to roll down- 
ward along the negative gradient towards the minimum. The computer 


method based on this idea would use the following correction vector: 


* 
oo Cee 
where 
* oy -1/2 
a EEE (A-10) 
Gj ~ 3a id? 


and d is the distance of the step along the gradient. Before the 
process is begun, a value for d must be found. An appropriate 

step distance is difficult to find, however, and a different value may 
be needed for each new iteration. Moreover,the correction vector, 

a , May Sometimes point nearly 90° away from the minimum (Figure 10, 
I) . In the method of scoring for parameters (39) , both difficulties 
are reduced by assuming that the M-surface, at point &, es 

quadratic in the form it would have if the peak model were linear. 

The first derivative of M (the gradient) and the assumed second 
derivative (the information matrix) are calculated and the minimum of 


the surface having these derivatives is used for the correction vector: 
xl * -11 
ied = =f g (A ) 

~u 


Often one side of the M-surface is steeper than the other. 
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Figure 10. Visualizations of minimization surfaces: 
I, effect of damping factor (schematic). The closed curves 


ace contours of constant value’of M. Line OA is in the direction 


of the negative gradient, ee ; ab point Os” OB ws the correction 


vector, ie ; and oOC is the damped correction vector, Ay . The 
dashed line is the locus of the damped correction vector as the damping 
factor, p , goes from zero at point O to SHnEIpLCy- ae DOr te n=. 


Il, effect of step along damped correction vector. The 


labels are as in I . The damped correction vector, line OC , can 


easily be extended to point D which is closer to the minimum than 


mort C . 
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If the initial estimate is on the less Steep side, the correction 
vector, At , May overshoot the minimum and arrive ata point even 
higher than the starting point. The reciprocal of a damping factor 
added to the diagonal elements of the information Matrix, he - 
gives a correction vector that is between the negative gradient and 


oe in direction and less than ne in magnitude (38) (Figure 10, I): 


-1 
A, = -({"+1/p U) g* (A-12) 


where iY is a unit matrix and p a positive damping factor. As 
the value of the damping factor approaches zero, the damped correction, 


Na , approaches the direction of the negative gradient, approaches 
P 


zero in magnitude, and the slope, dM/dp , approaches the square of 
| 2 


the length of the gradient, lg (40) . The best damping factor, 


the value that gives the lowest value of M , was found by a method 
similar to that used by Pitha and Jones (41). Starting with a trial 
value (p equal to 40 or the best value from the previous iteration), 


a second trial value, P. , was tried, which was the value that gave 
Ene sminimum Of a quadratic function of 1n(p+0.1) fitted to My (the 


initial value for which Po is equal to zero), My ;oand Ba M7 p 


Bip equal to zero. The third trial value was at the minimum of the 


quadratic fitted to Mo ; My , and M, . The fourth was fitted to the 


lowest three previous values of M In anomalous cases where the 


extrapolated value of p was less than zero, arbitrary values of p 


were tried that were fractions or multiples of the best value of p 


obtained so far. Further attempts after the second were stopped when 


Pp was greater than 10° or when the difference between the predicted 


Minimum of M and the actual value was less than 5% of the difference 
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between the actual value and a ae Eran M that was less than M 
could not be found after four attempts, the lowest value of Pp was 
divided by 10 until such a value could be found. Sometimes the 
optimal value of the damping factor was difficult to find. More than 
one minimum in the graph of M versus p occurred in some cases. 
Figure 10, I shows a minimum along the dashed line at polntasc and 


another between points C and oO. 


As the damping factor becomes small, the damped correction 
vector, bo , becomes small and the decrease in M also may become 
small. The rate of convergence was greatly imporved at little cost 
in computing time if a one-dimensional search was carried out for the 
lowest value of M along the line defined by the damped correction 
vector, Ap , When the angle between g and Sp was less than 50° 


(38) . 


=a. (A-13) 


where d is the length of the step. The optimal value of the 
stepping factor was found by a method similar to the method used to 
find the optimal value of the damping factor except that a quadratic 
function of the stepping factor itself was used for extrapolation and 


interpolation. Sometimes a step of 1000 times the length of the 


damped step could be made. Such a step could be the equivalent of 


: ‘ 7] e 
Many iterations using the damped step alone. The circumstances wher 


a step along A is valuable occurs when the M-surface has a steep 
E 


i ar 
face that flattens out sharply, as when Ky and Ke Se ae ae 


the critical region (Figure 10, II). 
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Minimization processes can be terminated in a number of ways. 
One of the more common methods is to stop when the greatest change of 
the values of all of the parameters per iteration is less than a 
specified constant. The basic criterion of a minimum, however, is 
that the gradient (9) is zero at the minimum and the rate of change 
of the parameters per iteration has no relation to the gradient. The 
best criterion for accuracy is to stop when the length of the gradient 
is less than a constant. An economic criterion, as well as an 
accuracy criterion for termination, based on the cost of computer time, 
is usually necessary. One such method is to stop when the rate of 
decrease of M per second of computation time or per iteration is less 
than a given constant. Premature termination may occur with this 
method when the rate of change of M per iteration is low ina 
particular iteration due perhaps to difficulty in finding the optimal 
damping factor. The influence of a few poor iterations on termination 
can be minimized by adding fractions of the rates of decrease of M 
per previous iterations to the rate for the current iteration. We use 
the following algorithm. Starting with Do equal to My ; DS Lor 


iteration mn is calculated from the formula 


where f is a positive constant less than one (0.1 was suitable). 


| i less 
Then, if lg*| was less than a constant, k, ; OD LE DM was 
\ 


than a constant, ky , the process was terminated. The value 0.0001 


was suitable for both constants. 
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Computer programs were written in Fortran IV and run on an 
IBM 360-67 computer with the Fortran H compiler; 96K of storage was 
required. The computation time for a single peak varied from a few 
seconds to nearly one minute; in general, the computation time was 
inversely proportional to the width of the peak and increased when 
the fit was poor. The computation times for the computer-generated 
peaks, for which a good fit could always be obtained, were usually 
about half the time required for a real peak with similar parameter 
values. The median computation times for normal heptane and normal 


metane were. _4.5. and..4.2 sec. 
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APPENDIX III 


COMPUTER PROGRAMS 


CONTENTS* 
Main Programs Page 


ii Estimation of. Peak Parameters, 1, 6.a, 7.a, 7.b, 
ee Sei eels) pe 18 fy De ca GW Sm an cell Ge ae oe Se ee ce Se Oe law cen 128 


MZ. ‘Retention Index. Measurement, .1, 6.a, 7.a, 7.b, 
Peele nebo, 7 ,. 18, 21 PRES Be a? Sat 4. Wee sek ee Oo PR ee Lot 


M3. Computer Generated Peaks, Lye OL age casy Tb eo, 


BEES LO, Le) (8 GS a! ZO GM ke OO cae oe ced eel co tee vee ae Gee 134 
Mew frower—Density Spectrum Estimation, » 12,.16,. 17 «. . « «. 139 
Cee OMA OVS is Og Ss Loe 4a cerlde oe set aw ce Ke Ge Se a oe Oe ee 141 
Subroutines 
ee: = ea Ano. Di peotCymeoypeles Gl AY ch wo Ge Sw ee 144 
x, INFORM ee AIR A Bee Tae, Ak CN CSS RR cs oto Re ee 149 
oo TNE SSS) ee ee et re ee eee he er fae ere ee 152 
PSE Doe Gey ein ee eh Tat ce i, SG Se se Sees oes 155 
Se  DSIMO Ra: ok es Ten Be SES COMET inl a beleer Gs) eu nee SemeOuaES one 158 
PEDO AMES io, 6 ee RP ee dk les. eet ero eee 160 
6b. TEST 
6c. PRIOR 
hast SET era eke Boop AS ee cc@) oD ae 5® 38 (b8 Lo ye¥ 'o® e85 OR Cem orem sees 162 
7b. CONVRT 
8. DPINIT er Oy ck ce eas ees 2a, 8 Le G6 he eee eee 163 
165 


9. FTPEAK (A,X) ° e . ° ° ° e e e ° ° ° ° ° ° e ° e ° ° ° ° 


Srey see 166 
MME EAS (AUN) i+) GE. ce ie Se Se ees 


“The numbers listed after each program or subroutine are the 


subroutines called within that program. 
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Library Subroutines 
Zoe CH201IA = provides date and 


22s CS003A,C - provides random 
distribution 


time of day 


number from 


uniform 


23. CSO19A,B,C - provides elapsed computation time 


Bane CSO1L2ZA —- matrix inversion 


Output Examples 
Poo remeM 7) $f. LO «Ee oe 


Pog ram Msp emt fF Geb file) 
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otra | 


Wee) 


180 


det 
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PROGRAM M1 
BAYESIAN ESTIMATION OF PEAK PARAMETERS 
SINGLE -PEAKS 
PARAMETER ESTIMATES ARE LISTED IN THE ORNER: AREA, 
RETENTION TIME,WIDTH, THEN THE THIRD TO SIXTH CUMULANTS., 
THE PARAMETER HAVE DIMENSIONS OF THE INPUT DATA. 
IMPLICIT REAL*8 (A-H,0-Z) 
INTEGER NAME(17),1P(5),NDATE(6) 
REAL*4& CGRAM 
COMMON/DATA/X(512),W(257),CGRAM( 2000) 
COMMON/ PARM/A(7),Y(7),V(7,7),PA(7),PV(7) 
DIMENSION DIF(512),DX(257) 
N=512 
CALL CH201A(NDATE) 
GAEL SSE? 


READ IN SET INFORMATION AND NAME, 
CARO Mia shea Pe?) lees) ,HLEVEL, FTES 
PPGLYCIX,PLI-OUTPUT PARAMETER 
=0 -MINIMUM OUTPUT 
=5 -=MAk IMUM OUTPUT (SEE BAYES) 
teGe) Chl je = STATISTECS OQUTPUT “PARAMETER 
=). =9NO OUTBUT 
=UPORge = GRAPH OR SERROR CURVE 
=3 - AS FOR 2 PLUS REAt PARTS OF TRANSFORMED 
DATA FANDGERROR CURVE (SEE BAYES ALSO) 
£205). Clipe a-ANALYSTS CASES -GSEE DERE 
=0 NO PRIOR’ INFORMATION 
=] PRIOR INFORMATION -SMALL PEAKS 
=2 PRIOR INFORMATION -LARGE PEAKS 
HLEVEL(F8.0)- PARAMETER FOR DPINIT 
FTEST(F10.0) - TERMINATION PARAMETER FOR BAYES 
END OF SET MUST BE FOLLOWED BY A BLANK CARD. 


READCSMIOD) IIR CL) PIP C2) ple Gs), HEEVEL, FTEST 
PEEFTESTsERsLJ00<06) ETEST=1. 00-05 
LFCHEEVELLLTS1.0D00) HLEVEL=5. 0000 


READ IN INVERSE OF POWER-DENSITY SPECTRUM (SEE PROGRAM M4) 


READ ES 105) CWO), J=1, 257) 
CALL) OCLPRT(257,W) 


READ IN CHROMATOGRAM. 
LAST CARD OF CHROMATOGRAM MUST BE FOLLOWED BY A BLANK CARD 


[FP LASTOCARD )RUEL. 
10 CALL DATARD(CGRAM,M, IERR,NORUN) 


READ IN PEAK Po cee 
CARD 1- NRUN,NOP,L1 , NAM 
NRUN(12)_ - CHROMATOGRAM NUMBER 
NOP(12) -PEAK NUMBER 
1P(4,5)(214) - SAMPLE NUMBERS BETWEEN 
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WHICH PEAK !S LOCATED 
NAME(17A4) =NAME OF PEAK 
CARD Q=| PRIOREYMEANS PAC1),1=1;7  C7FEI0.0) 
CARD 3- PRIOR STANDARD DEVIATIONS | 
Pil jaw leedle oie Go ceO. 0 ) 
CARDS 2 AND 3 ARE OMITTED WHEN 1P(3)=0 
THE LAST SET OF PEAK INFORMATION MUST BE FOLLOWED 
BY A BLANK CARD. 


70 READS, 102 Je NRUN, NPoeLP Ghd, LPC5), (NAME (1), P=1, 177) 
LEC WP CS) CEOs 0) 7 GO TO 10 
DROP)? Gls DJ ANOn VERR.EO.0) WRITE(6,999) 
WRITE(6,202) NORUN,NP,I1P(4),1P(5),NAME,NDATE 
GOe gd 0. 2-2 

PE were Piso). GT.0) READ(5,;104%) PA,PV 
GOP tO" 20 

PPLE CUERReGl.0),GO TO 21 
e205.) EO 10) "GO: 10. 24 
READ Co, LON) PA, PV 
Welle tbh, 20). PA, Py 

24 |Ph=IPC(4) 


OBTAIN INITIAL ESTIMATES 


CAGEADENECT CIP, HLEVEL) 
HaGhe Ch). La. 0) GOn TO, 20 
behPOs bP Ch) +1 


CONDITION THE DATA 
CALL CONVRT(L) 
OBTAIN ESTIMATES 


CALUBBAYES GIP, TEST, Cer SM,SEC, 01F) 
FEEPPO2).EO.0)°GO.10 20 
FeCl PG). bt.3) GO T0238 
DO 26 a =25,256 
[=N+2-d 
DX(J) =DSARTCX( J) **24+XC1)**2) 
26. CONTANUE 
DX(257)=NSQRT(X( 257) **2) 
omen iil FP PPEAKCA,%,01F) 
Pepe Cope 5) GO 10 0 
DOS 2 sda2 7257 
DIF(J)=DX(J)-DIF CJ) 
52, CONT LNUE 
DIF(1)=0.0D00 
DX(1)=0.0D00 


PRINT OUT REAL PART OF FOURIER TRAMSFORM OF CHROMATOGRAM 


CALL DCLPRT(257,0X) 
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C TABULATE REAL PART OF ERROR CURVE 


GCEEDEUCEERI C2o7, 0 Le) 


PRINT OUT GRAPH OF ERROR CURVE 


50 


30 


rigs 
102 


104 
105 


202 


204 
99:8 


a9 


CALL DIFFT(9,X) 

TEMP=10.0D000*A(3) 
[F(A(3).GT.10.0D00) TEMP=100.9D00 
IF(A(3).LT.2.5000) TEMP=39.0N00 
KS=!IDINT(A(2)-TEMP) 
KL=IDINTCA(2)+TEMP) 

Reeds. UT. 1) KS=1 

FEOKL IGT OL) KEEL 


L=KL-KS+1 

DO 30° J=1,L 
[=I PL+KS-1+d 
K=kS+d 


TEMP=CGRAM(1) 
QU SuENgMP RCO 


CONTINUE 

WRETE G6., 9.9.3-) 

Gime GRAPHCX, ULE )A, Vo Woks 200,29 
GOO. 2.0 
BORMAGGLY,.5.01..58...0 .F 10...0,) 
EORMATC212,214,1744) 
FORMAT(7F10.0) 

FORMAT(8D10.5) 

PORMALGCLIHO pteneRUN, wabnl 2s. BREAK 
FORMAT,CIX,.2 P7D15...5.) 

FORMAT(1HO) 

FORMAT(1H1) 

END 
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PROGRAM M2 


BAYESIAN ESTIMATION OF PEAK PARAMETERS 
MULTIPLE ISOLATED PEAKS FOR RETENTION INDEX MEASUREMENT 
PRINTS OUT SUMMARY OF RUN 


IMPLICIT REAL*8 (A-H,0-Z) 

REAL*4 CGRAM 
COMMON/DATA/X(512),W(257),CGRAM( 2000) 
COMMON, PARM/ACT) 7 YC7),V07,7),PAC7),PVC7) 


CLE S VON Det?) APCS, 7),-SOAC5, 77,105) 5.5D 15), ARCS) 


il POUR oT OS) 

VATEGER. PAMEXT) VPS 9ENPTCS):,NDATE CE) 
N=512 

CALL CH201A(NDATE) 

CALL SET 


PROGRAM PARAMETERS SAME AS FOR PROGRAM M1, EXCEPT REAL 
PAR OOF DATA AND ERROR? CURVE’ NOT PRINTED FOR] IP C2) 43 


Reese) srk a le C25 1 P03), HEEVEL Piesa 
Pewee vec.isls W000) HLEVEL=5.0000 
Peceiror. Lt. 1200-06) FIEST=l. 0D=03 

READS, 2001) (WCU) d=1,257) 

DAL DCUPRI(257,W) 


10 !0K=1 
CALL DATARND(CGRAM,M, LERR,NORUN ) 


READ IN PEAK INFORMATION 
PEAK NUMBER (NP) - 1 = AIR PEAK, 2 - SOLVENT, 
3 = NORMAL HYDROCARBON, & - COMPOUND WHOSE RETENTION 
INDEX 1S TO BE MEASURED, 5 - NORMAL HYDROCARBON 
IF 1P(5) FOR SOLVENT PEAK IS BLANK, 
SOLVENT NOT ANALYZED 


20 READGS5 1100.2) NRUN,NP,IP(4),1P(5),NAME 
HROIPCUSD.£O.0) GO TO 10 
HWAIGPE(I2EQ:0) GO TO 4&1 
WRITE(6,999) 
WRITE(6, 2000) NORUN,NP,1P(4),1P(05),NAME,NDATE 
GO? TO" 245 
22 PRAGHPAS)) (G10 ) REAND(5,1003) 
GO 1FO" bd 
2y GE GNERR.GT.0) GO TO 22 
IFCIP(3);.EQ.0) GO TO 26 
READ'CS 2008) PA, PV 
Were (6, 2002)..PA,PV 
26 CALL DETECTCIP,HLEVEL) 
reoppamlr.0) GO’ TO 40 
L=IP(5)-1P(4)+1 
NPTCNP)=L 
CALL CONVRT(L) 
IPL=IPCh) 
GAUL BAYES(IP,FTEST, P,FSM, SEC, DIF) 
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Pe OUPK ATE. 3) GO TQ 29 
CBELTETPEAKCA,X) 
CALL-DILFFT(9,X) 
KS=IDINTCAC2)-10.0000*A(3)) 
Sebo TCAC2) +10 20000*A(3)) 
bACKS ETI L) KS=1 
PECK. Gis. LY KLet 


L=KL-KS#1 

DO 28+ Jai, 4 
P=|P4+KS-Ll4+J 
K=KS+J 


TEMP=DBLE(CGRAM(1)) 
DIFCJ)=TEMP-X(K) 
28 CONTINUE 
GALEAGRAPHOX;DIF,A,Y,W,L,300,2) 
29 *ReTP CH) ~LT 20) GO-TO 40 
A(2)=A(2)+DFLOAT(IP4)-1.0000 
DO=S0VU=1,7 
AP(NP,J)=ACJ) 
SDACNP,J)=YC(uU) 
30 CONTINUE 
TCNP)=A(2) 
SDTC(NP)=V(2,2) 
FS(NP)=FSM 
PeGNP wet): GO" tO “20 


[IF ERROR OCCURS, SET SUMMARY DATA TO ZERO BUT CONTINUE TO 
ANREYZE OTHER SPEAKS 


GO'T0:50 


40 !0K=-1 
41 DO 42 J=1,7 
AP(NP,J)=0.9D00 
SDAC(NP,J)=0.0D00 
42 CONTINUE 
NPT(NP)=0 
SDT(NP)=0.0D00 
FS(NP)=0.0D00 
TCNP)=0.0000 
AR(NP)=0.0D00 
SDR(NP)=0.0D00 
GO TO 20 


CALCULATE SUMMARY DATA 


50 IFC IOK.GT.9) GO TO 60 
DRS 13 =17,5 
T(J) =0.9D00 
SOT(JU)=0.0D00 
AR(J)=0.0000 
SOR(J)=9.0000 
52 CONTINUE 
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RETI=0.0D00 133 


SDRET=0.0N00 
GO PEO 70 
Benen 
DO go uSt)s 
fFCNPT(J)-2E0.0) GO TO 62 
ACJ) =SDT(J) 
TVET CH= 10 
SD RGUNSDSOR TESHTG )+500) 
AGd) SAP Chat) ARGS, 1) 
SDR(U)=DSQRT(CSDACU,1)/AP(J,1))**24+(SPA(3,1)/ 
1 AP(3,1))**2)*AR(J) 
61 CONTINUE 
62 CONTINUE 
SDT(1)=0.0D000 
SDR(3)=0.0D00 
XN=DLOG1O(T(4)/T(3)) 
wreO0G10(T(5)/TC3)) 
RETI=1.0D02*XN/XD+7.09D02 
R1=XD/TC4) 
R2=XN/T(5) 
R3=(XN-XD)/T(3) 
SDRET=DSQRT(R1**2*AC UY )+R2**2*A(5)+R3¥*2*A(3) 
: +(R3+R1-R2)**2*A(1)) *1.,0D002/XD**2 


70 WRITE(6, 2003) NORUN,NDATE 

NP=k 
WRITE(6,2004) NP,RETI,SDRET 
WRITE(6, 2005) 
HA IME CON 2006 POOLE CAROL) 3Y, dere T em 1), ARC), 
1 Ne RChieteal os) 
WRITE(6, 2007) 
ipimr Gon 005 Leas WN, (SDAC IJ), J=1,7),SDECIO7 SOR Cl), 
1 ESC) lat,o? 
Go LT0b10 

999 FORMAT(1H1) 

1000 FORMAT(1X,311,F8.0,F10.0) 

1001 FORMAT(8D10.5) 

1002 FORMATC2Z12,214,17A4) 

1003 FORMAT(//) 

1004 FORMAT(7F10.0) 

POOCREORMATCIHO,', RUN *7102,' PEAK “pb2 2G oA) 


2002 FORMATC' PRIOR INFORMATION - MEANS AND STD. DEVS.' 
20021 /2(4X,1P3015.5, 0P4F10.6/)) 

2003 FORMAT('1DATA SUMMARY RUN NUMBER ',12,80X,6A4) 

2004 FORMAT('O RETENTION INDEX PEAK NO.',12,' =',F10.5 


20041 hOX, LSHOy DENEPSlI E1076) 

2005 FORMAT(! PEAK PARAMETERS! 74X, "NET RETENTION RELATIVE 
20051 AREA NO. POINTS!) 

2006 FORMAT('0',13,1P3D15.5,0P4F10.6,F14.4,F14.38,110) 

2007 FORMAT('OPARAMETER STANDARD DEVIATIONS BEX; 

70071 hcliM. OF SQUARED: DEV: ~) 

2008 FORMAT('0', 13,1P3D15.5,0P4F10.7,F14.5,F14.9,1PD14.4) 


END 


a es NY ; : tu) Seem 


i= 
wl 
a 


(hG2+CL E82) FADE 
| 


‘ 
* eke bs pe eee ; \‘\e ¥ 7 
4 PE ie J \ Pe ie fo be} BA 7 
aera Seo? ase Te, oe 
( 4! ‘ aca , ye. 2 re er 45 
. 14” £ \, \ ~ : 
} ‘ r i me ; 
opus 
> Sura 


CLEDTA (2 paeaaaes 
LGU a AeON +¥ é ref OROvd= IT 
. 4 hone 
RTA 

. LEP eeenx 

Ae tine t hat PA «hee 2a ae T ite 


i 
KV DoT rent het Be 2 ee ic Roe | 7 
de ; oe 


é : f ; 


° 4 i i a va " . bs iba ae OF 
ie ' | et ; ere at 0 | Naa 


. (ihe a eRe 
t ava 1 .oTe (TA ee mr 


(una (i a, St. a 


Z013,%2 1 Shy ms 
“SV TAP ir ce 


(A118 f9 
a aaa 


is, om ‘ 
Rial a  ' . ah 
_ Venn sTAOL a 


i ree wa 2. 


Re ee 


- 


CUO: Go) GaGa.) 


Go Qs. 


GaGa aqaan 


@@ © 


134 

PROGRAM M3 
ANALYZES COMPUTER GENERATED PEAKS ANDEN TO REAL BASE 
LINE NOISE. AS MANY TRIALS AS NESIRED MAY BE RUN ON THE 
SAME, SEGRION OR CBASE LINE, HE AR THIETC DALR PEAKST ARE 
GENERATED FROM EDGEWORTH'S SERIES EXTENDED TO THE TENTH 
CUMULANT.  BASEat INBeSEC TONS LiliMil T BDSM. +5 12PIP.OWNTS), 
AN INTERFERING PEAK WITH THE SAME SHAPE AND WIDTH AS THE 
PRANCTPAL PEAK. CAN. BE INTRODUCED... 
THE SIMULATED CHROMATOGRAM CAN BE MULTIPLIED BY WHITE 
MOUSE] TO SIMULATES Jd PERS 
ERROR, CURVE OF. LAST, TRIALe IN SERIES-4S PLOTTED 

IMPLI CI Ta REAL*8. (A-H;0-Z) 

INTEGER 1P(5),NDATE(6) 

REAL*& CGRAM,RET,DOWN,UP,RN(512),ALPHA 

DIMEN SIONSEXC5d 2oy ATA dd, DAC7T),DVC7),VMC7),VVC7), 

1 Pol NitOsbrde2edes X22 S27) 

COMMON PARM/Alvay, YGi oe C7, 7)7,PAC7),PVC7) 

COMMON/ DATA/X(512),W(257),CGRAM( 2000) 

N=512 

CALL SET 


READ IN GENERAL PROGRAM PARAMETERS 


RAs 10.01) Phe hee?) tPC3S.) ABLEVEL, BREST 
FECHEEVELS LT 31.0000) .HLEVEL=5. 0000 

LpUaeec tans e00-06) FTEST=1.00=03 

PEC het? oki. dil RG2)=1 

UP baie CT) 

IP2=1P(2) 


READ IN INVERSE-POWER-DENSITY SPECTRUM 


READ(5,1002) (W(J),J=1,257) 
CALLA DCLPRE(25a,W) 


READ IN BASE=-LINE NOISE 


5 CALL DATARD(CGRAM,M, 1ERR,NORUN) 
LECLERRAGT.0) STOP 1 


READ IN TRUE VALUES OF PARAMETERS. (TWO CARDS) 
RETENTION TIME, A(2), MAY BE BLANK 


10 READ(5,1003) (AT(C1),1=1,7) 
LFCAT(1). EOeg0e0D00) 760980 5 
READ(5,1003) (AT(I),1=8,11) 
eCup eem..0) GO a0 11 
READ(5,1003) (PACI), 
READ (5.2003) CPVC1), 

YP) CALL CH2O1ACNDATE) 
IPC1)=1P1 
LP G2) ee 
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NEXP = NUMBER OF TRIALS TO BE MADF 
INITIAL = INITIAL RANDOM NUMBER FOR RETENTION TIME 
GENERATOR, (POSITIVE INTEGER) 
[Ph, 1P5 = NUMBER OF DATA POINTS LEFT AND RIGHT OF 
RETENTION TIME TO BE USED BY ESTIMATOR 
TOTAL NUMBER OF POINTS USED = IP4+IP5+1 
SECMAX = IF TOTAL ELAPSED COMPUTUTATION TIME FOR THE 
EXPERIMENT EXCEEDS SECMAX, FURTHER TRIALS WILL 
NOT BE MADE, BUT THE SUMMARY DATA 
WILL BE PRINTED OUT 
DOWN, UP = RANGE OF RETENTION TIMES, 
ALPHA = PROPORTIONALITY FACTOR FOR MULTIPLICATIVE NOISE 
NO NOISE ADDED IF ALPHA=0.0D00. 
Z = INITIAL NUMBER FOR NORMAL-WHITE NOISE GENERATOR 
(MULTIPLICATIVE NOISE) 
AS1 = AREA OF INTERFERING PEAK RELATIVE 
TO PRINCIPAL PEAK 
AS2 = DIFFERENCE BETWEEN RETENTION TIMES OF INTERFERING 
AND PRINCIPAL PEAK IN UNITS OF 
WIDTH OF PRINCIPAL PEAK 


READ(S,1004) NEXP,INITAL, IP4,1P5,SECMAX, DOWN, UP, ALPHA, 


PECEXE. Fo. t SMER DER ana 


NEAP SG NEXP= 
VAR ND EG6;1 9.919") 
WR BP TEMEGt £2:010.0) CADPHALIZ,, AS1,AS2,NOATE, CATCID bel, 7), 
a NEXP, UNT TIAL, 1-P4, TPS (SECMAX (DGWN, UP 
WRVP ECG 2 005 2jCATCI), 138,11) 
ReCuR@GIr.GT.0) WRITECG, 2006) PA,PV 
GALL COS.C.03FANCiERTAL) 
TIME=90.0D00 
BOY 12701 =1517, 
DA(1)=0.0D00 
DV(1)=0.0D00 
VM(1)=0.0D00 
VVC 1)=0.0D00 
CONTINUE 
FM=0.0000 
FM2=0.0D00 
PEXP=1 


GENERATE WHITE NOISE IF REQUIRED. 


1£(Z.LT.1.0000) GO TO 20 
CALL RNG(N,RN,Z) 


DO 14 J=1,N 
RN( J) =ALPHA*RN (J) +1. 9 


14 CONTINUE 


GENERATE RETENTION TIME. 
20 CALL CSOO3C(DOWN, UP, RET, NRAND) 


AT(2)=DBLE(RET) 
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PA(2)=AT(2) 

IPCY)=1IDINTCAT(2))-1P4 

EUG eA LEW Dl a i es Bl gd bo 

IPCS)=IDINTCAT(2))+1P5 

PeCEPRS ECT. ND) se TRUS) =N 

ber Pt =P Ch) et 

oe Cle Gi ed) AND. CIEXP.GT<1)) WRITE( 6,999) 
GENERATE TRUE PEAK 


Caer PePEAK(CAT XSEX 
PEC AC 19ers ty OD00)*GO-40 10 
Ped exec, NEXP) CALL DELPRT(257,/EX) 
CALLE DT EET (OS) 
A(1)=0.0000 
A(2)=0.0D000 
LFECAST. GT. 0.0000) GO TO 210 
DO 205° J=aTSL 
P=1P(4)-1lt+d 
PINTCJ)=DBLE(CGRAMC I) ) 
205 CONTINUE 
GOrto: 225 


GENERATES INTERFERING PEAK- IF REQUIRED 


210 AC1)=AS1*AT(1) 
A(2)=AT(2)+AT(3)*AS2 
DQ 215 J=5,7 

ACJ) =AT(uJ) 
225 CON TUNUE 
CALTEETPEAKCA, PINT) 
CALE. DVEE TCS, PINT) 
HDO° 220° J=1, 4 
[=1P(4)-1+d 
PINT(J)=PINTC I) +DBLE(CGRAM(1)) 

220. CONTINUE 
NO2) =AC2)+1-1P CH) 

225 A C29 =AT(2)4+1-1P C4) 
WRITE(6,2001) PEXD PAT(2) -RETP IPOH) ARC) e ACL) Ae? 
[FCALRHATGTL 02091 GO TO 235 
DO 250ed=17t 

a1 PCP pede) 
TEMP=PINT(JU)+XC1) 
X(J)=TEMP 
EX(J)=TEMP 

230 CONTINUE 


CLOG ORL E 
235900 240 J=l, Lv 
P=7PCh) t0-1 


REMP=DBLE(RN(J)) 
TEMP=(PINT(J)+XC 1) ) *REMP 
X(J)=TEMP 

EX(J)=TEMP 
PINT(J)=PINT(J) *REMP 
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240 CONTINUE 


C 
; OBTAIN INITIAL ESTIMATES (MODIFIED FORM OF DETECT) 


Z50 CALLYDETECTCIPLHUEVEL) 
Ch PCS VETO. GOt TO) 20 
CALL CONVRT(L) 
Pe CVEKPR NE. NEXP)1GODT02 255 
PP Clijas 
1P(2)=3 
gop TAELIBAVYESCIPZFIEST,P,FSM,SEC,X2) 
POPC. LT.0) GO TO 32 
Ge S0rh=157 
DACT)=DACH)+ACI)-ATCI) 
DVCT)=OVC1)+C ACT) -ATC1)) ¥*2 
EXP=Y( 1 )**2 
VMCL)=VMC1)+EXP 
VVCT)SVVC TL) FEXP**2 
30 CONTINUE 
FM=FM+FSM 
FM2=FM2+FSM**2 
2 WME =TPMETSEC 
VAGwIiME, GT. SECMAX) 1G0) TO.50 
Pe Cle Cay. C130) GO 70-20 
peCraxrsGesNEXP) “GO TO: 50 
PPOCTeCLI.GT.0)" TEXPHLEXP41 
GO7TO! 20 


PLOT ERROR CURVE FOR LAST TRIAL 


CiGoG 


SOLCALL FIPEAK(CA, X) 
CALL DIFF T €9, XD 
KS=IDINT(A(2)-10.0D00*A(3)) 
KL=IDINT(CA(2)4+10.0D00*A(3)) 
PeCKStt i, to. ho=4 
IFCKL.GT.L) KLe=L 
L=KL-KS#+#1 
DO~28  J=1, Lb 
l=KS+J-1 
XCJ)=EXC1)=XC1) 
PINTCJ)=PINTCI) 
K={+1P(4)-1 
EX(J)=X(JU)-PINTCJ) 
28 CONTINUE 
CAML MGRAPHCEXEROPINT SA, Y,L, 200,99 
C 
: PRINT OUT SUMMARY OF EXPERIMENT 
PEC texP.1.2) GO: 10: A0 
EXP=DFLOATCIEXP) 
EXPT=DELOATCIEXP-1) 
DO S27 l=1,7 ene 
DACT)=DA 
Be PGnIGi EEL ASDA Cleo) /EXE! 
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VMCL)=VMC1)/EXP 
VVCL)=(VVCT)-EXP#YMC1) ##2)/EXPL- 
¥C1)=DSORT(DV(1)) 
VVC 1) =DSQRT(VV(1)) 
IF(DV(1).GT.1.00-20) DY(1)=VMCI)/DV(1) 
PACT )=DSQRT(VM(1)) 
IF(YC1).GT.1.00-20) AC1)=DAC1) *DSORTCEXP)/Y(1) 
52 CONTINUE 
FM=FM/EXP 
FM2=(FM2-EXP*FM**2)/EXP1 
1FCFM2.GT.1.0D-20) FMS=(FM-DFELOAT(IP(5)-IP(4)-7))/ 
1 DSQRT(EXP*FM2) 
TIME=TIME/ EXP 
WRITE(6,999) 
WRITE(6,2000) ALPHA,Z,AS1,AS2,NDATE, (AT(1),1=1,7), 
1 NEXP, INITIAL, [Pl, 1P5, SECMAX, DOWN, UP 
WRITE(6,2002) LEXP,FM,FM2,FMS, TIME 
WRITE(6, 2003) 
WRITE(6, 2004) (C1,DACI),YC1),DVC1), ACT), PACL), VMC1), 
1 vV(1)),1=1,7) 
GO TO 10 
999 FORMAT(1H1) 
1001 FORMAT(1X,311,F8.0,F10.0) 
1002 FORMAT(8D10.5) 
1003 FORMAT(7F10.0) 
1004 FORMAT(2110,215,3F5.0,E10.2,F10.0,£10.2,F5.0) 
2000 FORMATC(1HO, "EXPERIMENT PARAMETERS',1PE15.2,0PF15.0, 
20001 1P2D15.2,20X,6A4/ 
20002 4X,1P3D15.5,0P4KF10.6,13,18, 214, 3F6.0) 
2001 FORMAT(1X,'EXP NO.',1%,1PD22.5,0PF10.5,2110,1PN22.5, 
20011 0PF15.5) 
2002 FORMAT(1HO, 'N=',13,2X,'SUM OF SQUARES - MEAN =",1PN15. 


20021 5, VARPANCE =',015.5,' =' 015.5, 
20022 UMIAVERAGE TIMECSEC.) =',0PE10°1) 

2003 FORMAT('0',1T25, 'OBSERVED', 185, 'PREDICTED'/ 
20031 TOM DEX SIMEAN PDE VE Nee, "STD. IBV IPRaX, 


DOE O REET CIENCY’ -7%, USTUDENTS T', 15%, STD. DEV tak, 
70035 Vavneisnp. oF VAR. ') 

2004 FORMAT(1HO,13,1P4D17.5,7X,3D17.5) 

2005 FORMAT(4U9X,4F10.6) sans eeerene nen 

2006 FORMATC!" PRIOR VALUES'/4X,1 .5,0PHKF10. 

20061 ‘PRIOR STANDARN DEVIATIONS '/4X,1P3D15.5,0P4F10.6) 


END 
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PROGRAM M4 
ESTIMATES POWER-DENSITY SPECTRUM 
GIVES PUNCHED CARN OUTPUT FOR USE IN PEAK ESTIMATION 
PROGRAMS (WITH AND WITHOUT GROUPING CORRECTION). 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*4 CGRAM( 20900) 

DAMENSWONEX(C5 12) ,UX€512),WC512) 

N=512 

N2=256 

K=257 

PA=Oe 28513530/7179586/512.0000 

Z=1.0000/510.0D00 


READGANCONELSET Oba DATA.(512 DATA POINTS). 


h CALL DATARND(CGRAM,M, IERR,NRUN) 
PECIVERRSGTSON CSTOP 1 


Clune Rewhtited /2,1,~-1/ 2) 


X(1)=0.0D00 
L=N=1 


DORGUIU=2, 1. 
X(J) =NBLE(CGRAM(J)-0.5600*(CGRAM(J-1)+CGRAM(J+1))) 
6 CONTINUE 
X(N)=0.0D00 


FOURIER TRANSFORM 


CAREADEET (CS, xX) 
DO 8 J=2,N2 

br) Fah J 

Ns ie eae ee eee 
UX) = Khu) 
8 CONTINU 

XCK)=CXCK) **2) *Z 
UX(K)=X(K) 


SMOOTH WITH HANNING SPECTRAL WINDOW (1/4,1/2,1/4) 


VAR=X(2) 

W202) 30.75 *X02) 40. 25%XC3) 

DO 10 J=3,N2 
UAR=X (J) 
X(J)=0. 25D00*(VAR+X (J+ 1) )+0.5D00*UAR 
VAR=UAR : 

10 CONTINUE 
X(K)=0.5000+*(XCK)+UAR) 


CORRECT FOR FILTERING AND GROUPING 


DO 12 J=2,K 
HC=(1.0D00-DCOS(FA*DFLOAT(J-1))**2 


X( JU) =X CJ) /HC 
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UX(J)=UXCJ)/HC 
WCJ)=X(C JU) *(FA/(2.0D00*DSINCFA*DFLOAT(JU-1)/2.0000)) 
zk?) 


12 CONTINUE 


20 


ACL Y=XC2) 
UX(1)=UX(2) 
WC1)=W(2) 

Wel be (6777000 ) 
GAEL -DCLEPRTCK,.UZ) 
WRI TEG6, 2001) 
CALL OCLPRTCKS XxX) 
WRPTEC6,\2'0:0.2 ) 
GALL DELPRICK,W) 
DO 29 J=2,K 
X(J)=1.0000/X(J) 
W(J)=1.0D00/W(J) 
CONTINUE 
X(1)=0.00D00 
W(1)=0.0D00 


PUNCH OUT DUPLICATES OF INVERSE OF NONGROUPED 
AND GROUPED SPECTRA 


1000 
2000 
2001 
2002 


WRITEC7,1000) (X(J),J=1,K) 

WRITEC7,1000) (X(J),J=1,K) 

WRITE(7,1000) (WCJ),J=1,K) 

WRITE(7,1000) (W(J),J=1,K) 

GO TO 4 

FORMAT(8D10.5) 

FORMAT('IRAW POWER-DENSITY SPECTRUM’) 
FORMAT('QHANNED POWER-DENSITY SPECTRUM! ) 
FORMAT(C'OGROUPED AND HANNED POWER-DENSITY SPECTRUM! ) 


END 
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PROGRAM M5 
DETERMINES CORRELATIONS BETWEEN ERRORS OF PARAMETERS 
ESTIMATED BY PROGRAM M2, 


IMPLICIT REAL*8(A-H,0-Z) 
OMEN STON COACH @b6,7),B07),SCk, 7),CC757) 


READ IN PARAMETER VALUES 
ee CODE NUMBER OF PEAK 
See NORMAL HEPTANE (16 RUNS) 
= 2 = NORMAL OCTANE (16 RUNS) 
Seo pore TR IME THY LPENTANE C11) RUNS) 
See TOLUENE (5 RUNS) 
S PEAK AREA=DIVIDED BY WEIGHT OF COMPONENT. 
SVNETY RETENTION TIME. 


READ(5,1000)° ((CACIND, I,J 
READC5S,,1000) {XAQS, 150 Cd 
READCS, 1000D)7 GOACE 1,0) 7u 


=15 16) (ND=175.2) 
11) 


NORMALIZE AREAS 


Be Vesela) .11 
AREA=(ACT, 171) +A(2,1,1)+A03,1,1))*1. 00-05 
RO aNGiy (LP AREA 
AC 26121) =A(2,1,1)/AREA 
BC 5) 47a) =AC3, 1, LIPAREA 


205 CONTINUE 


WRER=(ACI, 1,1) 240(2,1,1) +404, 1,1)0041.00-05 
ACTON IO SACl, LOA PARE 
AC25 1,1) 2AC2, 1,1) / AREA 
Acts te 1 =A(h 1,2) 7 AREA 
210 CONTINUE 
13=1 
}4=16 
IND=0 
1 IND=IND+1 
Gomto C15) 15515778; 200) 77ND 
13 |4=11 
Gomu, 15 
Th 13=12 
14=16 
OBTAIN MEANS 
PDO 2S 1 S157 
B(1)=0.0D00 
SC IND, 1)=0.0009 
DO 20 J=l1,/7 
C(1,J)=0.0D00 
20 CONTINUE 
25 CONTINUE 


WRITE(6,2000) IND 
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DO 35 1=13,14 
NWRITECS, is IPOAGIND A do) od=1,7) 
DOmS0 eels 
B(J)= Bd) +ACIND, re 
30 CONTINUE 
Bo CONTINUE 
DO 4O J=1,7 
B(J)=B(JU)/DELOATC 14-1341) 
40 CONTINUE 
WRETEC6, 2002) 


OBTAIN COVARIANCES 


qanm 


DOjsss 1=15,14% 
)0750>4)=1,5.7 
ACIND, 1,0) =ACIND,1,J)-BCJ) 
DO 45 K=l,Jd 
CLARK) =CC, KD AGING Ip U)eACIND ETO 


ys! CONTINUE 


a VTEC, 2001) TP OACIND, tod), u=157) 
55 CONT! Hg COP oe aie ae 


DO~60)rFJ=1,77 
SC IND, J) =DSQRT(C(J, J) /DFLOAT(14-13)) 
60 CONTINUE 
WRITE(C6, 2003) 
DO 70¢b=13,/14 
DO 76S.edsL 7 
ACIND Jk, doRAGIND phad 2 ASGINDAd) 
65 CONTINUE 
WRLTEC6,2004) L,CACINOD,1,J),J=1,7) 
70 CONTINUE 
WRITE(6,2005) 
WRITE(6,2001) IND, (B(J),J=1,7) 
WRITE(6, 2006) 
WRITE(6,2001) IND, (SCIND,J),J=1,7) 
WRITE(6, 2007) 
C 
C CALCULATE CORRELATIONS AND CORRECT FOR BIAS 
C 


DO 80 1=1,7 


J=1,1 
of "CCL, JY=CCl_ JOS CSCIND, LY eSCUND WD ADELOAT CER=13)) 


C(1,J)=C(1,J)*(1.90004(1.0D09-C(1,J) *#*2)/ 
: (2. OD00*DFLOAT(I4-13-3))) 
75 CONTINUE 
WONT GS. 2008). betCUlsdiy vel, 
80 CONTINUE 
Coe Boel 
100 WR 6, 2006) 
ian: 2001) CIND, (SCIND,J),J=1,7),IND=1,4) 
WRITE(6, 2009) 
IND=0 
106 IND=IND+1 
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WOVE SMOG MOUS 10 5e 2100 a FN D 


MOAT iVareil 
‘We? 
jee 
bh=16 
C070) 110 
O02 1236 
bicead a 
G0 510-110 
0 senlne 2 
GO TO 110 
LOU =1 
Wah 
oe 12 
hive 6 
GOPTOn 110 
105) dda? 
CALCULATE CORRELATIONS OF PARAMETERS BETWEEN 
DIFFERENT PEAKS ON SAME CHROMATOGRAM 
110 DOTHRQ0R1=157 
61) =0. 0000 
DOME Sis sot h 
ROIS dS Te ACI2 J, 1 +B CT 
115 CONTINUE 
BCI)=BC1)/DELOATC( 14-13) 
BC1)=BC1)*(€1.0090+(1.0N09-BCI)**2)/ 
1 (2, 0D00*DFLOAT(I4-13-3))) 
120 CONTINUE 
MRUTER OS 2 OO)? E15 12 
GO TO 106 
290 STOP 
1000 FORMAT(7F10.0) 
2000 BORMATG'IDATA SET',13, NORMALIZED AREA") 
2001 BORMATGNO%, 14 ,F17.0, Dee ) 
2002 FORMAT('O DEVIATIONS FROM MEAN!) 
2003 FORMAT(C'O NORMALIZED DEVIATIONS') 
2004 FORMAT('O', Th, TELT. 4) 
2005 FORMAT('O MEANS!) 
2006 FORMAT(C'O STANDARD DEVIATIONS') 
2007 FORMAT('O CORRELATIONS BETWEEN PARAMETERS OF SAME‘, 
20071 "SPEAK! ) 
2008 ORMATC'O', 10X,15, 7F 15.4/) 


2009 FORMAT(' OCORRELATIONS BETWEEN PEAKS') 


2010 


FORMAT('O', 
END 
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SUBROUTINE BAYES(1P,FTEST,P,FSM, SEC, X2) 
TIMATES PEAK PARAMETERS 
P(1) IS OUTPUT PARAMETER 
Ore reer vor tad. ESTIMATES, ‘SUMMOEUSOUARES 
DAMPING FACTOR, COMPUTATION TIME, 
BUT NO HEADING 
=lr= ASP FORF OF SWI tHP HEADING 
Seer hoe POR lige PLUS ESTIMATES AT EACH [TERA ION 
=Se=emAS FOR 2, PLUS OUTPUT FROM DAMP AND STEP 
FTEST = TERMINATION FACTOR 
Pv="DAMPING “FAGTOR® CRETURNED) 
FSM = MINIMUM 'SUM OF SQUARES! (RETURNED) 
SEC *= ERAPSED* COMPUTATION TIME CSECONDS) 
Me = "WORK SPACE %257) 


PMP ELUCIT REALS *CA=H,0=-Z) 

REAL*4 CGRAM 
COMMON/DATA/X(€512),W(257),CGRAM( 2000) 
COMMONPPARM/ ACT YsY C75 VO7TF77, PACT), PVC7) 
DIMEN STON XD USPMC?7 77), RHS C7 ),XX€7),BBC7, 705 BROT); 
L DVAMCTO FAM 7), ALC7)>X% 2:01), C C702 CM G7) 
PMmeGer £PC5),PIVOTC7]) 

DATASNOL2/ PN27256/,K/257/,D0ROP/90. 1 D007 
NTIME=2500000 

CALL CSO1LSACNTIME) 

NIT=0 

IND=0 

beG@te Gi)..Gi.0) WRITECS, 1000) 

DO =2sJU=Z77N2 

}=N+2-J 
K2 @) =i Cd) e OX C0) ee 24X01) e42) 

CONTINUE 

X2(K)=WCK) *XCK) **2 

CALL INFORM(0,FSP,X2,C) 

FSM=FSP 

DESENT=FSP 

GO TO 5 


OBTAIN INFORMATION MATRIX (V) AND GRADIENT VECTOR (Y) 


mM 


IN 


CALL INFORM(1,FSP,X2,C) 
TRODUCE PRIOR INFORMATION, IF ANY, TO V AND ¥Y 


Cet eale Gs). G1.,0) scALL PRIOR 


FSO=FSP 
GL=0.0D00 
SCALE PARAMETER SPACE TO IMPROVE NUMERICAL PRECISION 
DONGE \-i77 
AL(1)=AC1) 
AMCI)=AC1) 
DIAC1)=1.0D00/(DSQRT(V(1,1))) 
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YC1)=Y(1)*DIACH) 
GL=GL+Y(1) **2 
6 CONTINUE 
GS=-2.0D00*GL 
DO 9 I=1,7 
DO 7 J=1,7 
| VOI, J)=VC1,J)*DIACL) *DIACU) 
CONTINUE 
CONTINUE 
CALL CSO19BCKTIME) 
SEC=2.6D-5*DFLOAT(NTIME-KTIME) 
P=0.0D00 
BECTEC1)LLTL2) “GO TO. 10 
WRITE(6, 999) 
WRITE(6,1001) NIT, (AC1),1=1,7),FSM,P,GL,SEC 


oo Om! 


CALCULATE DAMPING FACTOR 


POPCALL DAMPCIP, IND, P/F SP,GS,GAMMA) 
PECENDSGE: 100) “GQ 10°30 
Bai, OD00+1. OD00/P 
PEHDOr 1s t=177 
VC IIeZ 
13 CONTINUE 
CRBUEESOL2ZACV; DISPM, RHS,XX,.7,BB,RR, PIVOE) 


CALCULATE CORRECTION VECTOR AND NEW ESTIMATES 


CL=0.0D00 
GAMMA=0.0D00 
DOtT2Z2¢ b=157 
CORR=9.0N00 
DOn2Z0KVvsl77 
CORR=CORR+Y (J) *DISPMC I,J) 
20 CONTINUE 
GAMMA=GAMMA+ CORR*Y(C1) 
CL=CL+CORR**2 
CORR=CORR*DIACI) 
AC 1)=ALC1)+CORR 
C(1)=CORR 
22 CONTINUE 
GAMMA=GAMMA/ DSQRT(GL*CL) 
IF(DABS(1.0D00-GAMMA).GT.1.0D-05) GO TO 23 
GAMMA=0.0D00 
GO TO 24 
23 GAMMA=DARCOS(GAMMA) *57.295779 


TEST IF NEW ESTIMATES ARE REASONABLE 


Zin CALE LESH PRD 
BEG (EGS). ak Ts O}TGO FTO 760 


OBTAIN 'SUM OF SQUARES! FOR NEW ESTIMATES 
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CALL INFORM(2,FSP,X2,C) 
LE CESP.GE.FSM) GO TO. 20 
GAMAM=GAMMA 


FSM=FSP 

DO 26 1=1,7 
AMCI)=AC1) 
CM(1)=CC(1) 


26 CONTINUE 
REPEAT WITH DIFFERENT VALUE OF DAMPING FACTOR 


GO TOCLO 
30 NIT=NIT+1 


STEP ALONG CORRECTION VECTOR IF NECESSARY 


IF CGAMAM.LT.50.0D00) CALL STEP(1P,FSM,FS0O,AL,AM,CM, 
a X2,A) 
LEGEP CHD SLT .0) GO TO 60 


ENDSOF LTERATION! =*TEST*FORMCONVERGENCE 


32 IND=1 
FSP=FSM 
DON Suite l, 7 

ACT) =AMC1) 

5% CONTINUE 
DESENT=DROP*DESENT+FESN-FSM 
VWECDESENTZESN. GTZETEST },.GOITO )37 
WRITE(6,1004) 

GO TO 40 

57 WEGGECGTI<FTESTITGO TO 35 
WEVFEC6GL012) 

G0( tO. 40 

S5R1F GRIT. ET. 40) GO FenS6 
IPC 4) =-2 
WRITE(6,1005) 

GO TO 50 

36 VEC SECs LT..6050N00)71GOETO & 
IP(4)=-5 
WRITE(6,1006) 

BSLiL0CS0 


STALLS UbCS SECT LON 
[PG2y 1S STATUSTICS: OUTPUT PARAMETER 
=0 - PRINTS OUT DEGREES OF FREEDOM (NU), 
NORMALIZED CHI**2 (D) 
=1 - AS FOR 0, RETURNS DISPERSION MATRIX (V) 
=9 «AS FOR 1, PRINTS SCALED GRADIENTS, SQUARE OF 
SCALED GRADIENT LENGTH AND LENGTH, STANDARD DEVIATIONS, 
AND RELATIVE STANDARD DEVIATIONS. RETURNS 
STANDARD DEVIATIONS IN Y AND 
CORRELATION MATRIX IN V. 


ie 
h24 
fs} 3 
. (1) )MA 
‘ HO 
) O BTI 
OL 
(+T Uy 
, ) OAGAG 
/ ; 
, ra. tal bh oar 
L ee 
Arr C0,T (eA 
RY et 
: ~ = 
; Out BE 
23+929 |i 
aa se oe 
: ; pwRe t : de 
aur 
1at= fT} | eS OAtl= 
7 ee TRH a2 ay 


ha 01 09) 
26 OF Op CTeRTaII aa) A4 i 
a i (S9O1 a) aT IAW ies 


ae te Ot OT. on ; 
ar ot Of) POaVPD TIAL & 
i b--f0091 - a 
: Rioin se 7 
Z : | OT Op 
¢. GT OF {O60 +1. pa 


Ne on - 
Bata ae ih 
f sad a) aria 
‘ 2 cree 


notTaae 2012: TA 
AT aM TAS. TUSTUO SIITETTATS 2 


tUey MO aa 1°30 28% Lqysyrt THO 2THEAS ya 0 J 
| ka) Bs x0ha H25) JAAS The oa 
~ ntze An De +7) .0 AO] 2A i yay 
‘ST OATLO usa 42-87 Th I He. 3A - fe ; ry 
WATE JH prt Ry 1) ae a ts #2 [TARD a ror “ 


‘OLTAINGG CRATE GVETASER GHA 

ava ¥ Wa SUD TMYSO AROMAT: take 
ae ‘V UP RTARTA Oo pel pesriis are 
vie ‘gt ft ra. b.' 2a 


7 : : A 


kO 


4h 


45 


46 
48 


4g 


59 


60 


147 


=3 - AS FOR 2, AND PRINTS OUT CORRELATION MATRIY. 


Cur C2 yee Loy eee 10 56 
DO 44 f=1,7 
VCl,1)=1.0N00 
CONTINUE 
CAEL CSOLZACY FOTSPM/RHS PXX P77 BB RR PIL VOL) 
DRG ees il had ard g 
XXC1)=YC1) 
YC!) =DSQRTCDISPM(1,1)) 
CONTINUE 
DOP Ler = 1 77 
DO 46 J=1,7 
WETS ed PED SPM eS) CrCl) * ¥ Cd )) 
CONTINUE 
CONTINUE 
G5 =07 0000 
Come our = 1] 
XXC1)=XXC1)*YC1)*2.0D00 
VOD) = VC PD TACK) 
PECDABSCACI)).GT. 1.0045) RRCIJ=YCl)/DABSCAGCI)) 
GS=GS+XX(1)**2 
CONTINUE 
GAMMA=DSQRT(GS) 
CALI CSOLSBCKTIME) 
SEC=2.6D-05*DFLOAT(NTIME-KTIME) 
WP ECG, 7001) NIT CAC!) , 1=2,7)5FSM,P,GL, SEC 
be Gio, Gi 02—GO 10. 60 
PE Ghee. ci. 2)" GO 10 50 
NUSIPC5)=1PC4)<7 
DEV=1.9DO1 
IFCNU.GT.10) DEV=DSQRT(2.0000*FSM)- 
DSQRT(NFLOAT(NU*2-1)) 


1 
WRITE(6,1013) XX,GAMMA,GS 


WeUTE CG. 1007) NUC DEY; CY Cl) Ts15 7) 
WRITEC6 1008) (FTEST,RR 

BEGleG2) lot.) GO) LO 60 
WRITE(6,1009) 

MPITECG, 1010) (IVC, J), 321,70, 9=1, 2) 
CAUL SCs L9C 

RETURN 
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999 FORMAT(1H ) 

CORMAN CROUET XP ATT 13K, PAE 3X LOA BX 
POOL AS Bh MEAG U, SX, NAT, LOX, ESM! 16X,1 PO, 9X, TIME ') 
MOOT FORMAT (1X, 15,1P3D15.5,0PUF10.6,1°D15.6,D10, 2,012.2, 
10011 OPF7.2) 

1004 FORMAT(10X,'TERMINATED BY SLOW CONVERGENCE’) 

1005 FORMAT(10X,'TERMINATED BY ITERATION LIMIT") 

1006 FORMAT(10X,'TERMINATED BY TIME LIMITATION) 

1007 FORMAT(' ESTIMATED STANDARD DEVIATIONS', 60X,'NU=',15, 


10071 5 Mee DS 5 Ks 5015.15, OPS PLO. 5) 
1008 FORMAT(' RELATIVE STANDARD DEVIATIONS (PERCENT) ',66X, 
10081 Ped eoie= urd Os / OX p2 Por lower Deo 


1009 FORMAT(' CORRELATION MATRIX") 

1010 FORMAT(1H0O,20X77F10.4) 

Poll ECORMATCS9X, 1PD15.6,010.2) 

1012 FORMAT(10X,'TERMINATED BY SMALL GRADIENT’) 

Pits EORMAT( '* SCALED GRADIENTS '/&X, 1P3D15.2,4D10.2,2012.2) 
END 
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Chr yehe) Og ape i ae «wht + 
HO TATE MELE: Ste MR Ee ae 2 TANn0t'R 

W! Xia, Sea SAARMATS ONAN Tee Bn 
‘ VE WHI ete EF we. erat. ae us a 
2NOY PALE) ARCHES av TASS, "ST AMAG 
ve. SL ees ata, T2273") Soe 
| CCN TAN A RELE iag TAN 
| LOOT IS CREE ‘GH EDTAY ! 
Ch APD, 0). CRORE AORTA Ha 
(1742 PARA. AR OFZ ARIMORT! MOLITA 
At Par XO Pe TT OAR. Basie) + ae 


MAO 


COG Gio GoiGa Geo Ga GaGa 


SUBROUTINE INFORMCITP, PHI,X2,C) 


IP 1S OPTION PARAMETER 
=0 - RETURNS INFORMATION MATRIX AND GRADIENT VECTOR 
AND SUM OF SQUARED DEVIATIONS 
=] - RETURNS INFORMATION MATRIX AND GRADIENT ONLY 
=2 - RETURNS SUM OF SQUARED DEVIATIONS ONLY 
X2 = WEIGHTED ABSOLUTE SQUARE OF DATA, USED IN 
CALCULATION OF 'SUM OF SQUARES' WHEN CONTRIBUTION FROM 
FITTED PEAK IS NEGLIGIBLE (REAL PART OF NORMALIZED 
Be Ne Loe Leos, THAN Exe te 20). ) 


IMPLICIT REAL*8 (A-H,0-Z) 
REAL*4& CGRAM 
COMMON/ DATA/X(512),W( 257), CGRAM( 2000) 
COMMON A PARM/ACT VN C7 jo VC7,7),PACT),PVC7) 
DENSON Er C7) OFC] )X2¢1),C(1) 
PRA N) S127, 1279567 ,M/7/, TWOR1/6.28318530/179586/ 
Pere eo s:2). GO 160. 100. 
DO 4& [=1,M 
Y(1)=0.0D00 
DO 2 J=1,M 
V(1,J)=0.0D00 
CONTINUE 
CONTINUE 
FA=TWOPI/DFLOAT(N) 
FM=FA*(A(2)-1.0D000) 
FS=FA*A(3) 
PECIPUNESLDPPHT=0. 0000 
DO 20 J=2,N2 
S=DFLOAT(J-1)*FS 
$2=S*S 
F=((-A(7)*S2/30.0N00+A(5))*S2/12.0D00-1.0N00)* 
1 $2/ 2. 9D00 
PEE GT .-20.0000). GO-TO 7 
IFCIP.NE.1) PHI=PHI+X2(J) 
GO TO 18 
7 PECRSLE  PeOD01y"GO TO 6 
WRETE CH; 2000) YCACTE) (1 Us1,7) 
GO TO 40 
8 ALPHA=DEXP(F) 
$3=S82*5S 
SA=DFLOAT(JU-1)*FA 
=N+2-J 
ARG=(=AC6)*52/20. 0D00#AC4) )*S3/6.0D00-DFLOAT(J-1) 


OM — ho 


CF=NDCOS(ARG) 
SF=DSINCARG) 
DPR=ALPHA*CF 
DP1 =ALPHA*SF 
10 PR=X(J)-AC1)*DPR 
P1=X(1)-A(1)*DPI 
Tar eMN ees PH| =PH1+U*(PR*PR+PI*PI) 
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sat oF-8 cS. 03.41 359 | 
Wes ed iat 8 ad 
(0,G=€1)7 = . 

eal 5.00 
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"tag 
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ANE TAK X= 8S 
* dweu 


14 


13 


A oy) 


eG Paik Qe GOs LOMTR 
EF(2)=0.0000 
0} tansy, \ 


EFC3)=((-A(7)*S2/20.0D00+A(5))*S2/6.0000-1.0N00)* 


S*SA 
OF(3)=(-A(6)*S2/12.0D00+A(4))*S2*SA/2.0D00 
EF(4)=0.0D00 
OF(4)=S3/6.0D0 
EF(S5)=S2*S2/2h.0D00 
OF(5)=0.0D00 
EF(6)=0.0D00 
OF(6)=-S3*S2/120.0D00 
EF(7)=-S3*S3/720.0D00 
OF(7)=0.0D00 
Y¥(1)=Y(1)+U*(PR*DPR+PI*DPI ) 

U=U*ALPHA 
T=U*ALPHA 
Meter yey C1.) )4+T 
DO FT? (K22 Mm 
YCK)=Y(K)+U*(PR*(EECK) *CF-OF (CK) *SF) 
+P] *(OF(K)*CE+EE (KY) *SF)) 
VOK,K)=VCK, K) +T#CEFCK) *EE CK) +90F(%) *OF (K)) 
CONTINUE 
DO 14 K=h,M 
V(3,K)=V(3,K)+T*(EF(3)*EF (CK) +0F(3)*0F(K)) 
CONTINUE 
V(1,3)=V(1,3)+T*EF (3) 
V(2,3)=V(2,3)+T*90F(2)*0F (3) 
Ho or kK=opie 2 
1K=K-1 
ee K)=V(1,K)+T*EF CR) 
1K) =V02,1K)+T*0F(2)*0F (IK) 
sour QUEUE 


ONCON TINUE 


Ze 
26 


28 


Phe ls GOS TO! 26 
PH! =2.0D00*PHI+X2(N2+1) 
PEChEP. 6Q. 2) GO 10 40 
V(4,6)=-V(5,5)*4.0000/5.0D00 
V(5,7)=-V(6,6)*5.0000/6.0000 
¥(1)=2.0D00*Y(1) 
ale is) =) ODUUSNCL,.L) 
FA=2.0D00*A(1) 
FM=FA*A(1) 
DO 34 1=2,M 
YC 1) =FA*YCI) 
VG1p I )=FASV CLS ID 
VG WO =V CP) 
V(1,1)=FM*VC1,1) 
K=!+1 
1eCK, COLOMF 1) 2G0 cl) 40 
DO 32 J=K,M 
V(t,u)=FM*VC I,J) 
V(J, 1) =VC1,J) 


Sif ot on ts $ Ae a 

. j ‘ a($}4a 

“~~ Ait oh i 

i 27H 4 HO10 OSA SEN IAS citi FAFpy 
«(Ce TACPARH LINGO OIA ys( 290) 
gph Ot tet & bat 

MOO NE 30), 

GHa0 PS cer eahaa” 

Oat ; Betz h30.. 

.G=( 48, : 

OGED , a cm og 

h OHO, OO Peek aent tage eo 
=e no0e. be(0)90 - 
| bye BG or cGr e ; i | Watt Ph a i. 

; Apa JAe Uet. o 


/ ROMIARUAT 


tT a 
A+ oo 


pli 
“a0 MEHR apa 
‘BUKIT HOS a 
BE Ot SRe a O3, qa 
(So S40) OX THA OGRO. So 1N9 | 
fd 61 ‘an KS.p3 layy4 mi 
Peet ea -4 0a Ve pully BY 
NPTG BAHL, 2+ (dob Weekt 2 a 
Ob ¥eOOG0. Lp ) 
in crt) Vantn.. Seeks vo 
\ a | GhIA* O00, iat | 


Anny see Te — 
fan aa sia. tispeates Hy 


Ve (fa ii 
a> id uw) : _ ' iain 
7 ) | - : ars 


f r 
> q 
ad 


CIC Gow GD. 


Be CONTINUE 
ate CON LL NUE 
4KO RETURN 


DETERMINE IF &'TH AND 6'TH CUMULANTS ARE 
Pier CAL REGION 


100 T=A(5)**2-1.6N00*A(7) 
Pelee UU) GO) PU. .b 
SA=1.0ND00/(TWOPI*A(3)) 
LA COABSGAL 7 ine d.00=10) .GO TO =b05 
GAs Lewin 10) GO 10.6 
S=DS0ORUCL2Z. 0000/AC5 ).)*SA 
Peo, Gil .,0 0 B00). GG 10 6 
GO. -f/0),.12.0 

105 T=DSQRT(T) 
S=(A(5)-T)*15.0D00/A(7) 
S$2=(A(5)+T)*15.0D00/A(7) 
Geo mee UP CU. AND. od, Le. 0. 0000) -GO°TO.6 
BEGG Gd. OOOO CAND 25 2..67.0.0000) "GO Ton tio 
S=DMAX1(S,S2) 
GOREG TTS 

110 S=NMINT(S,S2) 

115. S=DSQRT(S)*SA 
Fees. Gil.w 0000 GO" 10. 6 

120 .S=A(5)**2/1.6D00 
SOD sch) ACCT) 
AC7)=S 
GO TO 6 ) 

1000 FORMAT(4X,1P3N15.5,0PHF10.6) 
END 


15s 


er YA Weoens Tele Alse” 

CV TA onge,, e | 
OT Of CORRS ie eS CRA, 08 
OT GD COUGG se ith, ei OM a0 


WHS BAGO GO 


C 
C 
C 


Onn 


iG G 


SUBROUTINE DAMP( IP, IND,P,FS,GS,GAMMA) 


CALCULATES DAMPING FACTOR (P) FOR NEXT TRIAL 
UPC UTPUT PARAMETER, JIPCL), TS EQUALATO Ss ,PPRINTS OUT 
INFORMATION ON EACH TRIAL (SEE STATEMENT NUMBER 700) 


10 


IND = TRIAL NUMBER 
poe SUMLOF SQUARES FOR P 
GS = § 

GAMMA = ANGLE BETWEEN SCALED CORRECTION VECTOR ANN 


EOPE OF “ESGPR) WAT P=0 
NEGATIVE GRADIENT (DEGREES) 


PMOULCTTVREALSs —CA=H,0=-Z) 

DEMENS VOM ACS 7.8. G5) 7X C3). ,F C3) 

PNEEGERS RPGS } 

DATA AN/S/ /UP/ Fo 0D01/7, D0WN/ 1.0D-01/,H/ L.0D= 017 ; 
his 2. 5,0 25a 092094 ON b/ 


“HECIND JG. 0) "Go Sto Lo 


PL=40.0D00 
GO TO 100 

IFCIND.EQ.1) GO TO 100 

FSP=FS-FO 

[F(FSP.GT.0.0D00) PBAD=P*0.95D00 
IFCFSP.GE.FSM) GO TO 20 

FSM=FSP 

PM=P 


P0eGO TOSG1007200,,.300,400, 700,700), 1ND 


FIRST TRIAL - PREVIOUS VALUE 


100 -FO=FS 


FSM=0.0D00 
IF(PL.GT.1.90D07) PL=PL*DOWN 


PM=PL 

P=PL 
PBAD=1.0020 
F(1)=0.0000 
X(1)=HL 
IND=2 

GO TO 300 


SECOND TRIAL - MINIMUM OF QUADRATIC 


200 REZ YSRSP 


=DLOG(P+H) 
0 es WRITE(6,1000) GAMMA,FSP,P,X(2) 
Bide CED SAC 1) -H¥GS*(X(2)-X(1)) 1 / XC) XC TD) ee? 
B(2)=H*GS-2.0D00*B(3)*HL 
B(1)=F(1)-B(2)*HL-B(3) *HL**2 
IND=3 
GO TO 500 


THIRD: TRIAL 


00 TEGP SGT. b/0N05y' GO’ TO "700 
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507 


519) 


LSS 
F(3)=FSP 
X(3)=XP 
WeGhe GL) SGh 32) WRITECG, 1001) -BC1),B(2)5803),FP, 1c, 
L GAMMA, FSP,P,XP 
IND=4 
GO TO 459 


URTH- TRIAL 


Pace. Gist.0005) GO T0700 
WeGieCl)<GT. 2) WRITEC6, 1001) BC1),8(2), 813), ER ic, 
1 GAMMA RESP SPO XP 
GAMMA = ANGLE BETWEEN SCALED CORRECTION VECTOR AND 
NEGATIVE GRADIENT (DEGREES) 
Boer (1) 
K=1 
PoC?) tec e) GO TO 405 
FP=F(2) 
K=2 
CFGS9eLEsER) nGO,» T0420 


O 460 J=1,N 
ACJ) =1.00D00 
K=J+N 
ACK) =X(J) 
K=K+N 
ACK) =X(JU)**2 

CONTINUE 


LCULATE CONSTANTS OF QUADRATIC 


Chia BS iAMOCA, BN, 1) 
ALYZE CASES AND TAKE APPROPRIATE ACTION 


IFCIC.GT.0) GO TO 700 
FP=0.0D90 

Pauets)-G).0.0000) GO. 10 510 
[F(B(2).G7.0.0).GO.70.505 
1C=3 

P=PM*UP 

XP=DLOG(P+H) 

GO TO 300. . 
IF(B(1).LE.0.0D09) GO TO 507 
Nees, 

€0. 10 503 

IC=h 

GO TO 515 
XP=-B(2)*0.5N00/8(3) 
1F(XP.GT.HL) GO TO 520 
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1C=5 

515 P=P*DOWN 
XP=DLOG(P+tH) 
GO TO 300 

520 FP=B(1)+B(2)*XP+B(3)*XPx*2 
hE GEPOLE. RSM) “GO TQ °525 
1C=10 
GO TO 700 

eZ amir (XP PRES20 70000) LG0S5TO1530 
Ic=7 
GO 704503 

530 P=DEXP(XP)=-H 
1C=0 
PeGP,Gi.),00=08) GO TO 800 
P=1--00-08 
1C=6 


PRINT OUT IEORMATION 
B(1-3) = CONSTANTS OF QUADRATIC 
FP = PREDICTED MINIMUM OF QUADRATIC 
IC = ANALYSIS CASE CODE - SEE PROGRAM ITSELF 
FSP = DIFFERENCE BETWEEN WEIGHTED 'SUM OF SQUARES' 
ATiPs0 AO AT THE CURRENT VALUE Obl? 
P = CURRENT VALUE OF DAMPING FACTOR 
XP = LN(P+0.1) 
700 PRGPEUPYSGT.2) WRITE(6, 1001) B(1),802),B¢3)5P, IC, 
1 GAMMA, ESP,P, XP 
1FCFSM.LT.0.0D00) GO TO 710 
P=P*DOWN 
IND=6 
feeP NE. 1.00-10) GO TO 710 
GO TO 900 
710 P=PM 
PL=PM 
IND=100 
Go To 900 
800 IF(P.LT.PBAD) GO TO 990 
P=PM*2,0D00 
IC=8 
[E(P.GE.PBAD) P=0.5D00*(PBAD+PM) 
arent Uppase6 A2p1 0.2) 
C95 POLS 6; 
vie Bor ae dwdy io 3015.2, 015.6, 13,0PF10.2q0epdsesy 2nd 06oD 


END 


-“— 


7 : 
, 
ty 
~ Pm | 
ay a 
\ 
an 0 J 
qi! 
e 4 ; ‘ 
qt fiAsM 
fit ep t 
\ C8, Oto Ph 
“ j - —" 
eva c 


| (2-98, 
TRCN =F 
B1eY J ARAS "ORs 
AATATIG ~.983* 
PAA $=4 TA 
AN THARS tum myx 


-— 
ooa0 Tu: MeV 


6 r-tin, Lad. ADA . 


39 


iubieresaa 2 

DCI 0TH: 
OT Of 
ate 
Re 1785. 


Sac i 
? 


MOAT TUO TH 


(1,049 MS qx) ? 


~ RAL BR 


WOO» q=9 
7 PS ae Sie 


BOR, HT. 00. 
mqeq © 

, aed 

} éoreant. 
0k OT Dd. 
(MAS9.TI, ay4t 

) aods. Seiad. | 
aie ‘geor <4 7 & 

9 cconaee Seve = ee 
MAUTIA 008 

OLA XOvyTA Any 0 
GCSE aT aAOF 


QOD OQ Oo Qa @ Gare 


qan 


qOanann 


CA 
CO 
1F 
PR 


IN 


10 


CA 
PS 
AT 
2 


15 
16 


Ryd 


155 
SUBROUTINE STEPCIP,FSM,FSL,AL,AM, CM,GX2,A) 


LCULATES OPTIMAL VALUE OF STEP (S) ALONG DAMPED 
RRECTION VECTOR AND RETURNS NEW VALUES OF PARAMETERS 
OUTPUT PARAMETER, “1PC1T),OBQUALUS 3, INFORMATION IS 

INTED OUT AT EACH TRIAL (SEE STATEMENT NUMBER 400) 


FSM = ''SUM OF SQUARES" AT S=#1 

FSL = "SUM OF SQUARES! AT S=0 

AE = PARAMETER. VALUES AT S=0 

AM = PARAMETER VALUES AT S=1 

CM = CORRECTION VECTOR Ay S=1 

Gre e=a WEIGHTED ABS WARE 0 OF DATA 
A = PARAMETER VALUES” AT OPTIMAL §S 


(MRT ITT REALS3 CA-H,0=Z) 

OTMENSHON ALC) CMC/),GX201) (P13), 2% C36 ¥507 O89) AM CI OV AG?) 
MITEGER?) dP C5) 

DATARN/ 7/6/57 ,UP7 1. 0001/, DOWN/ 1200-017 


POMALIEZES SV TRY Seen? 


SBAD=1.0D12 

1C=0 

XMAX=1.9D04 
P(1)=0.0D00 
X(1)=0.0D00 

PM=FSM-FSL 

P(2)=PM 

X(2)=1.0000 

SM=1.0000 

S=2.0D00 

IND=0 

IND=IND+1 
XMAX=XMAX*1.0D02 

LE ( SHiePSeSBAD): GO TO 12 
S=SM*2.0D00 

1C=7 

1F(S.GE.SBAD) $=0.5000*(SBAD+SM) 


LCULATE NEW VALUES OF PARAMETERS AND "SUM OF SQUARES! 
IS THE DIFFERENCE BETWEEN THE 'SUM OF SQUARES 


S=0 AND AT THE CURRENT VALUE OF S. 
DO 15 J=1,N 
AGHJEAL G4 S¥OM( J) 
CONTINUE 
GAUL CTESTONP) 


PRIENRG@inieLT.0) GO TO 405 
CALL INFORM (2,FS,GX2,Q) 
IS 
eee ces 0D00) SBAND=S*0.95D00 
feyec cE.FSM) GO 10: 30 
FSM=FS 
PM=PS 
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SM=S 
DO 20 J=1,7 


AM(J)=ACJ) 
20 CONTINUE 


20 GOTO (100, 200,200,200,400), IND 


FIRST TRIAL, S = 2. P€1=3) AND X(1-3) ARE STORED VALUES OF 
ATTEMPED S AND THE |RESULTING PS. ey 


HOO ME CIPC1),.GT.2) WRITE(6,1000) BS jis 


P(3)=PS 
X(3)=S 
GO TO 300 
SUBSEQUENT TRIALS 
Comm ee@rra.)).G1.2) WRITECG, 1001) VOL) 2Y( 202703), POL ICL PSS 
mecOnose eG (PS=-PQ)/PS).LE.1.00-03) GO TO A405 
REPMACEMSSGIVING GREATEST PS WITH CURRENT VALUE OF <S 


208 PQ=P(1) 
K=] 
UE CP C2) VEE TPO) MEO) TOMZE0 
PQ=P(2) 
K=2 
20) LEC Gy we. BO eGo TTO S2E5 
K=3 
25 BCR) = PS 
X(K)=S 


300 DO 310 J=1,M 
Q(J)=1.9N00 
K=J+M 
Q(K)=XCJ) 
K=K+M 
Q(K)=X(J)**2 
S10 CONTINUE 


CALCULATE CONSTANTS OF QUADRATIC 
CALL DSIMQ(Q,Y,M, IC) 

ANALYZE CASES AND TAKE APPROPRIATE ACTION 
PEGG Gre )4G0 7000520 
PQ=0.0000 
IE(DA8S(Y(3)).GT.1.0D-20) GO TO 330 
IC=5 
PCY (2).GT.0.0000), GO 10° 320 

315 $=SM*UP 


Boer 10 
320 Boe DOHN 


BRE 


get rae 0.003 008.008} 
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GOLTON’TO 

S=~0. 5D00#Y C2) /Y C3) 

BE CY C3) .GT.. 050060) GO TO °350 
BECS NG PSSM) “GOrTO 30 

1C=2 

CO*TO RSIS 

1C=3 

GO TO '320 

PECYC2) CLP-OCOROO) GO"TO S60 
1C=h 

GOrt0\320 

PE OCOMET FXMAX) (GOTO 370 

1C=6 

S=XMAX 

GO«TrO 10 
PH=VETIOSECY (2) +S 47035 )) 
GOTO 410 


NT OUT INFORMATION 

Cros)" = CONSTANTS OF QUADRATIC 

Q = PREDICTED MINIMUM OF QUADRATIC 

ANREYSIES CASE CODE (SEE PROGRAM ITSEEF) 
DIFFERENCE BETWEEN SUM OF “SQUARES AT S=0 
ANGOEAT CURRENT VALUE OFS 

= CURRENT VALUE OF STEP FACTOR 


S 


fe CRR ER) UGTAD) WREITECG POOR) YC1L),Y(20,YO3) 72 IC PSzs 
RETURN 

FORMAT CSSX%,1PD15.6,010.2) 

PORMAN C19 Xs 123020.7,015.6,13,015.6,D10.2) 

END 
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SUBROUTINE DSIMQ(A,B,N,KS) 


SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS A*X=B 
DOUBLES PREEDSI@N FORM OF SIMO IN LBM SCIENTIFIC 
SUBROUTINE PACKAGE VERSION ITI 
N - NUMBER OF EQUATIONS AND VARIABLES 
Pa Gh BY eNYMATRAX OF COEFFICIENTS. 
DESTROYEDeIN COMPUTATION. 
Ba=e@VECTORTOP ORIGINAL ‘CONSTANTS. 
GEPEACED UR Y=SOLUT HONA=ax% 
OP = NORMAL SOLUTION 
1 - EQUATIONS ARE SINGULAR, SOLUTION NOT FOUND. 


KS 


MMPEPCLHT REAL ¥8 4 A-H, 0-Z) 
DIMENSION AC1),BC1) 
TOL=9.0D00 
KS=0 
JJ=-N 
NO 65 J=l1,N 
JY=J+1 
JJ=JJtNel 
BIGA=9.9N00 
1THdd-d 
DQ 30 JsJ,N 
PJ=1T+] 
IFCDABSCRIGA)-DABS(ACIJ))) 29,30,30 
20 BIGA=AC IJ) 
IMAX=I 
omc Ort LNUE 
IF(DABS(BIGA)-TOL) 35,35,49 
gor os 1 
GO T0)5100 
4O [1lsJ+N*(J-2) 
1T=1MAX-d 
DO 590 K=eJ,N 
P1=lt1+N 
P2eUbe T 
SAVE=AC 11) 
ACI1)=AC12) 
A(12)=SAVE 
50 ACI1)=AC11)/81GA 
SAVE =8( I MAX ) 
BCIMAX)=BCJ) 
B(J)=SAVE/BIGA 
PEGI): 55, 70,95 
55 10S=N*(J-1) 
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79 


80 
100 


DO 65 IX=JY,N 

IXJ=10S+1X 

IT=J-1X 

NO 60 JX=JY,N 
IXUX=N*(JX-1)+1X 
JIX=1XUX 41 T 

ACIXISX) =ACIXSX)-CACT XJ) *AC JX) ) 
BCIX=BCLX)=-C BCU) *ACT Xd) ) 
NY=N-1 

IT=N*N 

DO 80 J=1,NY 

[A=1T=d 

IB=Ne-d 

IC=N 

DO 80 K=l,J 
BCIB)=BCIB)-ACIA)*BC IC) 
LA=1A-N 

IC=1C-1 

RETURN 

END 
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SUBROUTINE DETECTC(IP,HLEVEL) 
PROVIDES INITIAL VALUES OF PARAMETERS FOR RAYES 


REMOVE STATEMENTS NUMBERED 1990 TO 1004 
WHEN USED WITH PROGRAM M3. 


cia Gp 


IMPEVCIT  REAL*8) (A-H,0=Z) 

REAL*#4 CGRAM 

COMMON /DATA/X(0512),(0257),CGRAM( 2000) 
COMMON/ PARM/A(7),Y(7),V(07,7),PA(7),PV(7) 
INTEGER? 1 ROS) 

N=512 

M=| PCS). TPGh) +1 


C 
C TRANSFER DATA 
C 
VOCGe CES=TP CH) 
1001 DO 1 J=1,M 
1002 b=L1l-1+J 
1003 X(J) =DBLE(CGRAM(C 1) ) 
1004 CONTINUE 
L1=M+l1 
DO*>2\J3J=Lk1,N 


X(J)=0.0000 
2 BOUT ARY 
OC 


18! 
GO. F 00,200,300),KP 


GASERO TEND, 2PREOR’ . PNEQRMATLON 
USES APPROXIMATE AREA, RETENTION TIME, AND WIDTH PROVIDED 
BY BRIN OTSISETS CUMULANTS TO ZERO 


PILPG 


£00 CALBYOPENIECIP,HLEVEL) 
(aG eur) sk T. 0): GO-10) 51 
NO 11 J=4,7 
A(J)=0.0000 
11 CONTINUE 
GO TO 490 


C 
C CASE 1 PRIOR INFORMATION- SMALL PEAKS 
C USES PRIOR MEANS AS INITIAL ESTIMATES 
C 
200:1D0 21 “J=177 
A(J)=PACJ) 
21 CONTINUE 


GOsEO 25:2 
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AAAS toi cassie Daiichi: 7 
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USES APPROXIMATE AREA, RETEN 
PRIOR MEANS FOR THE CUMULANTS 


300 


So 
oe 


BD 


OR 


IHFORMATION- 


CACLIDPTNETALPZHLEVEL) 
CEGLEXSOQGEAMO) GO TO 51 
DOMNBIN dab, 7 

ACJ) =PACJ) 


CONTINUE 


DO“55/UN=1,7 
PER PViGUN.GT.050000) PVC U)=1. 0000/7 PVC) )s*2 


CONTINUE 


ENTRY@TEST GLP) 


SETS5 PP.C4) 


400 


42 


LA CAC LOR 
Pe CAG2):. 
FECA 2) 
DE CACS). 


REL ACS). 
GOntor51 


48) 
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LARGE -PEAKS 
TION TIME, AND WIATH AND 


-1 IF PARAMETER VALUES ARE UNREASONABLE 


~9,.0000) 
«0. 0D00)) 
~DFLOAT(M) ) 
sO. 000) 


G90 TO 4&2 
GO TO eh2 
GO TO 42 
GO TO 42 


~ (DELOAT(M)/6.0000)) GO: 10°62 


WRITE€69208) (I,AC1),1=1,7) 


IPC kh )=-1 
GO) toast 


ENTRY PRIOR 


INTRODUCES PRIOR INFORMATION 


59 
51 


201 FORMAT(1HO, '*****ERROR SURQOUTINE TEST-UNPEASONABLE!, 


20a1 


DO 50 J=#1,7 
Y¥( JU) =¥Cd)- CAC) -PACI) ) * PVC) 
V(J,J)=V(d,d) +PVCd 


CONTINUE 
RETUR? 


LeGSTIMATES' /'0",7C' 


END 
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SUBROUTINE SET 


C 
C CALCULATES FILTER AND GROUPING CORRECTION 
Spee yGROUPRTNGE CORRECTION INOT  REQUPREDOGFIETER CORRECTION 
C ONLY, FOR USE WITH PROGRAM M3) REPLACE STATEMENT 
6 NUMBER 1000 BY: 
: HC(J)=1.0D00/(1.9D900-NCOS(2.0900*VAR) ) 
C 
EMPLACIT: REAL*8 (CA-H,0-Z) 
REAL*4 CGRAM 
COMMON/DATA/X(512),W0257),CGRAM( 2000) 
DIMENSION HC(257) 
DRAM £5 L272 5 oy, K/257/,P1/35.191592655509795/ 
FA=P1/512.0D00 
HC(1)=0.9D00 
DOMZeU=25:K 
VAR=FA*DFLOAT(J-1) 
1000 HC(J)=VAR/ (2. 0D00*DSINCVAR) **3 ) 
2 CONTINUE 
GO FO8s0 
C 
ENTR¥eGONVRTCL) 
C 
C FILTERS DATA, FOURIER TRANSFORMS, AND APPLIES CORRECTION 
C 
M=L-1 
VAR=X(1) 
X(1)=0.0090 
DO & J=2,M 
FA=X(J) 
X(J)=X(U)-0. 5D00*(VAR+X(Jt1)) 
VAR=FA 
& CONTINUE 


X(L)=0.00D00 
CALLSDEFTC9,X) 
DO 6 J=2,N2 
P=N+2-d 
X( J) =X( J) *HCCd) 
X(1)=XC1)*HC CU) 
6 CONTINUE 
X(K)=X(K) *HC(K) 
10 RETURN 
END 
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SUBROUTINE DPINITCIP,HLEVEL) 


Bere S ANTTITAL ESPIMATES OF ACT) ,A(2),AC3) 
LEVEEP EIS HEIGHT AT WHICH INTEGRATION IS TO BEGIN 
Ler ankon OCCURS, SPPCRITIS2SET)TO -1 


1 


PMPLICIT REAL*8 (A-H,0-Z) 
REAL*4& CGRAM 
COMMON /DATA/X(€512),W(257),CGRAM( 2900) 
COMMOWAPARIMT A GMIAY(7),VC7,7), PACT), PVC7) 
INTEGER @hPCS) 
N=(P GSS hRC) +I 
SOL ds1,3 
ACJ) =0.0N00 
CONTINUE 


FIND MAXIMUM 


F=X(1) 

M=1 

UG T2ud=2,N 
hE CXGd) “ERek? GOUNO 9 
F=X(J) 
M=J 


9 CONTINUE 


2 


CONTINUE 


FIND END OF LEADING AND TAILING TAIL 


IN 


PeCGM ot.) OR. (MGT. (N=$9))9° G0" T0838 


B2= CXCMF1 4X (M42 +X OM43 4X CMFh )+XOM45 ) #0, 29000 
R1l= CXC M=$1)4+XCM-2) +X (M=-3 4X CMH +X CM=-5 ))*0, 2N00 
K=M+1 

L=N-6 

DOF 3° d=K,L 


B2=B2+0.2000*(X(J+5)-X(CJ)) 
LE COX (00-82), LESHLEVEL) GO TO4 
CONTINUE 
GOTO 8 
K2=Jd+3 
L=M-7 
DO: 5) JHl Zl 
K=M-J 
BL=B1+0.2090*(X(K=-5)-X(K)) 
her (COM Bl). LEaHUEVELI GO Te 
CONTINUE 
BOOS 
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CALCULATE LINEAR BASE LINE UNDER PEAK AND PARAMETERS 


6 


8 
10 


102 FORMAT(1HO,!*x*** ERROR SUBROUTINE DOPINIT - 


OZ 


K1l=K-3 
SEOPES(B2=81)7 BE DOAT (K2=K1 ) 
eK 5 
DOT 7ed=K,U 
P=X Cu) =SLOPE*DELOAT (J=K1)=B1 
AC1)=AC1)+P 
D=DFLOAT(J-M) 
A(2)=A(2)+P*D 
A(3)=A(3)+P*D**2 
CONTINUE 
A(2)=A(2)/A(1) 
P=A(3)/A(1)-AC2) **2 
bE GPE Gt. 00-20) AC3)=DSORTCP) 
A(2)=A(2)+DFLOAT(M) 
GO TO 10 
WRITE CG, 102) 
PPC) =-1 
RETURN 


"PEAK WIDER THAN DATA') 
END 
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SUBROUTINE FTPEAK(A,X) 


GENERATES FOURIER TRANSFORM OF A PEAK IN THE FORM OF 
EDGEWORTH'S SERIES USING UP TO THE SIXTH CUMULANT. 

m =o SEVEN PARAMETERS 

Ae SRPBAK MCS ITZ OPOTNTS) 


IMPLICIT REAL*8(A-H,0-Z) 
DIMENSION Seis A(7) 
DATA PN fh Sei2' fh) N2/256/, RAZST/ 
FA= 6 .283185307179586/512.9000 
FS=FA*A(3) 
FM=FA*(A(2)-1.0000) 
X(1)=A(1) 
DO LTO Rd =2,N2 
}=N+2-J 
S=DEROAT(C J=1)*FS 
S2=S*S 
P=0CsAG7 I ¥S2 30 SOD00+FAC5))*S2/12. 000041, 0D00) 
i: #52 2. 0D00 
PERE SG h £2 20'-0 O00 G0 “TO: 92 


X(J)=0.0D00 


L0'15"¢ 596 0D00 


1.0000) GO TO 20 
ALPHA=A(1)*DEXP(F) 
ARG=(-A(6)*S2/20.0D004ACK) )*S*S2/6.9000 

1 -NFLOAT(J-1)*FM 
X(J) =ALPHA*NCOS(ARG) 
X¥(1)=ALPHA*DSINCARG) 
6 CONTINUE 
10 CONTINUE 
09.0D00 
a ae soaaieweDnsccn Uc 2)).GT.0.0000) WRITE(6,10900) 
RETURN 
20 WRITE” (06,1001) 
XC) 2] .000000 
1000 Levee deeae WARNING = SUBROUTINE FTPEAK - ', 
19001 'PEAK 1S ALIASED') 
1001 FORMAT ( ' **x*** ERROR - SUBROUTINE FTPEAK - -, 
10011 "BAD PARAMETERS') 
END . 
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SUBROUTINE FTPEAKCA, X,Y) 


GENERATES FOURIER TRANSFORM OF A PEAK IN THE FORM OF 
EVGEWORTH'S SERVES USING UP TO THE SIXTH CUMULANT. 
FOR USE WITH PROGRAM M1 
W=REAL PART OF e®PEAK FOR OUTPUT PURPOSES 
(257) ORD MORE! POINTS) 


IMPLICIT REAL*8(A-H,0-Z) 
DIMENSTONW ACL), X(1),YC1) 
MAA NG Saw /ZN2/256/,8/257/ 
PAO 2 051855071795 86/512.0000 
FS=FA*A(3) 
FM=FA*(A(2)-1.9000) 
XC1)=A(1) 
Y¥(1)=A(1) 
DOT 1] J=27N2 
P=N+2-J 
S=DFLOAT(J-1)*FS 
S2=S*S 
PaaS omo2/ 50. 0000tAC5) )*52/ 12. 0000-1. 9900) 
1 RO2/ 20000 
FeCPs Ga. -20.0000) GO TO 2 
X(J)=0.0D00 
X(1)=9.9D00 
Y(J)=0.9D00 
GO TO 6 
2 PE CEMGT. 051010009! G0. T0720 
ALPHA=A(1)*DEXP(F) 
ARG=(-A(6)*S2/20. 0D00+A(H) )*S*S2/6.9D09 
el -DFLOAT(J-1)*FM 
X( J) =ALPHA*DCOS(ARG) 
X(1)=ALPHA*DSINCARG) 
Y (J) =ALPHA 
6 CONTINUE 
10 CONTINUE 
X(K)=0.0D00 
Y(K)=0.0D090 
TEX DABS (xtND INEDABSIOUGKER) 1.67: 0.00000 WRITE(6, 1000) 
GOMTOT 25 
DOAWRETE COBSOOT) 
X(1) = 0.0D00 


25 RETURN 
1000 FORMAT(' **x** WARNING - SUBROUTINE FTPEAK - ', 


199901 IPEAK IPS ALIBSED > 
1001 FORMAT ( ' ***** ERROR - SUBROUTINE FTPEAK - LP 
LO0d1 "BAD PARAMETERS' ) 


END 
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SUBROUTINE FEPEAK(A,X) 


ERATES FOURIER TRANSFORM OF A PEAK IN THE FORM OF 
EWORTH'S SERIES USING UP TO THE TENTH CUMULANT. 
USE WITH PROGRAM M3 


IMPLICIT REAL*8(A-H,0-Z) 
DIMENSION AC1),X(C1) 
DSR AON 52 Att? 256 £ Ak / 2577 
Res = One? OOO. 4 795867 512.0000 
FS=FA*A(3) 
FM=FA*(A(2)-1.0D00) 
X(1)=A(1) 
HOR10 eJ=275N2 
1=N+2-J 
S=DFLOAT(J-1)*FS 
S$2=S*S 
B= _O0C0C-AC11)*S2/99.O0D00+A(9) )*S2/56.0D00-AC7)) 
eS 27m WOOFTA C5) Y*S2/12.0000—1.0000)*S2/2. 00 
PAGES e42070000) 1GO TO 2 
X(J)=0.0D00 
X(1)=0.9D00 
GO> TONG 
HPCE SGT 0 2ND00)/FGO.2TO 20 
ALPHA=A(1)*DEXPCF) 
ARG= (((-A(10)*S2/72.9D00+AC8))*S2/42.0N00-A(6)) 
*S2/20. 9N00+ACK))*S*S2/6, ANDI-DELOAT(J-1) *FM 
X¥(J) =ALPHA*NCOS(CARG) 
XC 1) =ALPHA*DSINCARG) 
CONTINUE 
CONT PNUE 
X(K)=0.0D00 
LECDABSCXCN2)Y4DABS(X(K-2)).GT. 9.0009) WRITECE,-L00D) 
GO TG425 
WRITE SCG AL00)) 
KCL) Fa 20 F000 
Soran a xeeee WARNING - SUBROUTINE FEPEAK - ', 
Noe we YS ALIASED') 
FORMAT ( ' ***** ERROR - SUBROUTINE FEPEAK - iy 
"BAD PARAMETERS' ) 


END 


cat 


SO: MAND BHT 1, nied ae on 2 
THK TES HTMST. vt sai ii | 


ouen) age sas 
| Cie Wor 14 MIC 

j ‘ppiluaies dea eke “ Tht ov: 
enon, sPe\nneets OEE, eR 


000. aatsL, 
re | (4) 
ah oe 
. ' basen 
aRett- vITAGLAE 


CET YA-HOMO, AeA TOKE Peete “oe\tes tl fikeaa ee = 
A RAG OH COM one SINGS sR) A OO. OF Wee ee a 
: OT ha Coote. ee oe 

‘ OGG .08CUIK :. 


o 
ator a ns) 


et OF 04 _(ouanyaT9. a it 


CAdAKIGeC 
CCAIA-CanD,¢ (Xoae heietg’, cp \gzacat need By san _ 
ited t=1,) + EU ARPA ELE VAENOOD 20 te exh, 
DRA ERI 
CPAMIGTEO* A 
7 


COMOr AVATISH (OORE 0, TH OCR RRC kam 


«* = MARIS serrated > gen nian”. 
IA-21VRAS > 


.' = wAae32 IY TORE + eo wéonw, * > 
CSERITINARAT. 


Oo AeA e ere a 


oO 


o™71IVD®?” 


168 
SUBROUTINE DFFTCLOG2N, xX) 


FOURIER TRANSFORM OF REAL DATA IN DOUBLE PRECISION 
X(N=2**LOG2N POINTS) IS REPLACED BY ITS COMPACT STORED 
FOURTER TRANSFORM. THE REAL PART FROM FREQUENCY F=0 TO 
Peet SeslORED AS XO))) PROM) J=? TOCN/ 241). THE -TMAGTNARY 
PART FROM F=l1/N TO,1/2-1/N 1S STORED FROM J=N RBACKWARNS 
MOes=N/2+2. SINCE IX 1S REAL, THE IMAGINARY PART OF 175 
ERAN SFORMY TS ®ZERO AT°FSO"AND F=1/2 CIN UMITS OF THE 
SAMPLING INTERVAL). 


M =f ROM p SRAT " A ey tpt 
U. OF Ae, EDMON on Ri STEN Ge Goo oNNbE at es neat Gee 
UMPIETO!T REAL*8™C€A-H, 0-Z) 
DIMENSION X(1) 
INTEGER J,JN,K,KN,N,N OVER 2 
Pl=3.141592653589793 
N=2**(LOG2N-1) 
CALL DRS80(LOG2N,X) 
CALL DRSBO(LOG2N-1,X(1)) 
CALL DRSBO(LOG2N-1,X(N+1)) 
CALDSDOMR TET (LOG2N=1,X(1),X(N+1)) 
N OVER 2=N/2+1 
DO 100 J=2,N OVER 2 
K=N+2-J 
JN=JU4N 
KN=K+N 
XR=(X(J)+X(K))*0.5099 
XP=CXCIN)-XCKN) )*0. 5000 
YR=(XCIN) +X (KN) )*0.5000 
Y!=(X(K)-X(J))*9.5000 
ARG=PI*NFLOAT(JU-1)/DFLOAT(N) 
C=NCOS(ARG) 
S=DSINCARG) 
T=YR*C+Y1*S 
YI=Y1#C-YR*S 
YR=T 
X(J)=XR+YR 
X(K)=XR-YR 
XCKN)=X1+Y1 
XC UN) =YI-XI 


100 CONTINUE 


KO=XCIVEXANSL) 
¥R=EX C1) EXC N+1) 
X(1)=XR 
X(N+1)=YR 
RETURN 

END 
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INVERSE FOURIER TRANSFORM 
INVERSE OF DFFT, APPLICATION OF DEFT AND DIFET 
hNeTHE aORIGENAESDATAs 

MODIFIED FROM PROGRAM OBTAINED FROM DEPT. 
OPrA., (EDMONTON, WRITTEN BY G. SANDE. 


BOOUENCE RESULES 


a8 


10 


100 


SUBROUTINE DIFFT(LOG2N, X) 


IMPLICIT REAL#*8& (A-!,0-Z) 
DIMENSION X(1) 
P1=3.141592653589793 
N=2**(LOG2N-1) 
N OVER 2=N/2+41 
N2=2«N 
FAC=1.9D00/DFLOAT(N2) 
DO 10 J=2,N 
l=N2+2<J 
X(J) =FAC#X (J) 
X(1)=-FAC#X(1) 
CONTINUE 
DO 100 J=2,N OVER2 
K=N+2-J 
JN=J4N 
KN=K+N 
AR=X( J) 4+X(K) 
AL=X(KN)-X(JN) 
BR=X(J)-X(K) 
Bl=X( JUN) +X (KN) 
AQRG=PIL*DFLOAT(J-1)/DFELOAT(N) 
C=NCOS(ARG) 
S=DSIN (ARG) 
T=BR*C+BI*S 
BI1=BIl*C-BR*S 
BR=T 
X(J)=AR-BI 
¥(K)=AR+B1 
X( UN) =BR+AI 
X(KN)=BR-Al 
CONTINUE 
AR=X(1) 4X (N41) 
BR=X¥(1)-X(N+1) 
X(1)=FAC#AR 
X(N41)=FAC*BR 
CALL DMRETC(LOG2N-1,X(1),X(N+1)) 
CALL DRSBO(LOG2N-1,X(1)) 
CALL DRSBOCLOG2N-1,X(N#1)) 
CALL DRSBO(LOG2N,X) 
RETURN 
END 
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SUBROUTINE D MReAFTCLOG2N,X,Y) 


MIXED RADIX ONE DIMENSIONAL FOURIER TRANSFORM 
Pe DOUBLES PRECISION 


MODIFIED FROM PROGRAM OBTAINED FROM NEPT. OF PHYSICS, 


Sis 


Or he EDMONTON, WRITTEN BY G. SANDE. 

MMP Calls (REAIAa CAH, 0-2) 

DUMENS LON, XCldy YC1) 

ries ppt 2 3 NMG 

Pero iGr Cig), C5,.b0, 41,12, 135,R0,R1,R2,R3751,52,S3; 


INTEGER A,B,C,N,E,F,G,H,1,d,%,L,4,N,8S,CS,0S,ES,FS,GS, 
eae NS KG eo MS ALS BL, CL,DL,EL,FL,GL ML, IL Su, RL, 
Walt, Mi, SCLBDy, UCTS) 


Sa eNO GBS55(2 )), CCS,SC3)),(0S,SC9)) (ES, S050, 
Gr. Co» ECS, Say Daas, SC3)),01S,809)), Gis, 5010), 
KS pease) ep kn 85) GE) BOM, 051013) CAL, UCL), CBU, UC20)0; 
Slory bt yi pede Utes) ), CFL, UC6) 2, (GL,UC7)), 
SG eWicpns Aly UU GO yd Jily Ute), CKL,UCI1)), 

2 kee UC 120), CML Clay) 


N=2**LOG2N 

Pa auoGers LE.1) GO TO 500 
DO 400 K=2,LOG2N,2 
M=2**(LOG2N-K) 
Mu = *M 
DO 390 J=1,M 
ARG=6.283185307179586*DFLOAT(U-1)/DFLOAT (M4 ) 
C1=DCOS(ARG) 

SL=DSINCARG) 
C2=C1*C1-S1*S1 
S2=C1*S1+C1*S1 
C3=C2*C1-S2*S1 
S3=C2*S1+S2*Cl 


avi We 
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: : Ay 
h 


4. padi Ten Lith 7. he 
ae “ETT AN) wie ne tine wf Ne tte at Bi 
mba Sil CEL SteL: 
(2)2.23). OSC CElAPA) KS ie eRe eee 

of PP LY COE ee) YEN TE EM OC BIE pen ay 

he) ) I) (CELIO URS b ODP 2a Oe 2 ee 

Su DR) ee hd Ie vlan 

/ \ ry 1 Ree ae a ie | 

| (tel 


ww ) 


coy 10.30¢0e ge 


beatae A ae 


YOSo: 
aw 


, 
= 
= 


[ 
Mi 
Se Sit ve 


=" 


eae 
ait i, 


gy 
"| abe + aes a 


a Lae 


ed 


| 
2% 
a 
coy ; 
ed Aa 7 _ ve ral yey 
| at Lie 7 
- ie _ 


DO 200 1=Mh,N,Mh 
JO=]+J3-Mh 
J1=J04M 
J2=J14M 
J3=J32+M 
RO= X50) #x(J2) 
R1=X(J0)-X(J2) 
1N=Y(J0)+Y(J2) 
Pt =y¥ GI0)=YCU2) 
R2=X(J1)+X(J3) 
R3=X(J1)-X(J3) 
P2=Y¥CJ1)#Y(d3) 
13=Y(CJ1)-Y(d3) 
X(J0)=RO+R2 
VUJOJ=10F12 
Peeueed. 1) GO TO 100 
ACJ2) =(R1I+83 )*C14+011=R3)*S1 
yCo2y = Gli-R3S)*Cl-GCRIF13)%*S1 
XCJIT)=CRO=R2)*C2+C10-12)*S2 
¥CJ1)=(10-12)*C2-(RO-R2)*S2 
Pevoye ORL = PS )FCSHCl T+RS )#S3 
1359 =( 11 6R3)*65-C€R1I=13)*5S3 
GO™ T0200 
100 CONTINUE 
ACUS)V=AR1L+13 
Y¥(J2)=11-R3 
X(J1)=RO-R2 
YCJT)=10-12 
X(J3)=R1-1 3 
GJ so. Lb EFRS 
200 CONTINUE 
300 CONTINUE 
4O0 CONTINUE 
500 CONTINUE 
Pe CULOG IN EO. LOGIN/2s2) GO T0790 
00,5600 P=1,N,2 
O=xXits) +x P41) 
R1P=XC lp) XCi41) 
Posy Gl)+y Cit) 
ey DY. Ct+1) 
X( 1) =R0 
Y(1)=10 
XC 1+1)=R1 
YC t+1)=<11 
600 CONTINUE 
700 CONTINUE 
MS=N/2 
ML=N 
DO 800 K=2,12 
J=14-K 
S(J)=1 
Jy=SCJ+1) 
eae AT 21a Std = St tl oye2 
800 CONTINUE 
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AL= 
JJ= 


NN 
DO 
NO 
DO 
DO 
NO 
DO 
DO 
DO 
NO 
DO 
NO 
DO 


JJ=JJt1 
Ute CU Ev) GOnal, 00 
T=X (Jd) 


BS 
0 
900 
900 
990 
909 
900 
909 
9090 
200 
9090 
900 
9090 
900 
3999 


A=1,AL 

B=A,BL,BS 
C2ReCL CS 
D=C,DL,DS 
ED,EL,ES 
Pee eieers 
CEEUGibGs 
H=G,HL,HS 
PSH 1, LS 
Jeb wts0S 
K=J,KL, KS 
L=K,LL,LS 
M=L,ML,MS 


X( Jd) =X(M) 
X(M)=T 


T=Y¥ (dd) 


Y¥(JJd)=Y¥(M) 
Y(M)=T 
CONTINUE 
RETURN 


END 
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SUBROUTINE DPD RS BO(CLOG2N,X) 


REVERSE SUBSCRIPT BIT ORDER 


MODIFIED FROM PROGRAM OBTAINED FROM DEPT. .OF PHYSICS, 


OF A., EDMONTON, WRITTEN BY G. SANDE, 


REAL*38 T,X 
DIMENSION X(1) 


INTEGER JJ 


VATEGEREAL BRT CHOMERR, Ghani w, KALEM, NABSACS, DSHESTES GS, 


Pt PMO NO, ee eel OE g ek, DuZeis PEG HEA Laie kk, 
~EU,ME,NL, SCD), UC14) 


SRUPVATENCECRGS, SCZ) ),(CS,S(3)90,COS7SCH)O CES 7S 5D), 
memento), (452507 )), CHS, 5(3)) -C IS, S09) 05 dS, 5010))7 
meno, oURTLDORCLS, S(12)), (MS, S15) ), (NS, SCs), 


sCAL, UTD) eOBE,UC2)), 
ceour Uae). (COL UCK ID) CEL, UC5)), CEL, UCGQ9, 
sno, CGS). 01 UPUCs)) CoE, UCI0)), CRE, UCIT) 


GU NCI, 
), 
CCU 2) VACHE (159 )CNL, UCIN)) 


NS=2**(LOG2N-1) 

NL=2*NS 

DOPLOOSKS2 215 

J=15-K 

U(J)=SCJ+t1) 

S(J)=1 

LEC SGUt2).GT.1) SCJ) =SCJF1)/2 
CONTINUE 


AL=BS 

JJ=0 

DO 200 A=1,AL 

NO 200 B=A,BL,BS 
DO 200.C=B,CL,CS 
DO 200 D=C,DL,DS 
DO 200 (ESD7ELEES 
D0-200"F=E)FL,FS 
DO 200 G=F,GL,GS 
DO 200 H=G,HL,HS 
DO 200° 4H,18;,1S 
NO 200 J=l,JL,JS 
DO 200 K=J,KL,KS 
DO 2001L=K,LE,LS 


DO 200 M=L)ML,MS 

DO 200 N=M,NL,NS 
JJ=JJ+1 

fee (Jd CE PNTEGO 1708200 
T= 

XC JSPVEXCH) 

X(N) =T 

CONTINUE 

RETURN 

END 


, E, : a 
i : 7 ; fr wn 
cx nombre “2 a wiTuoRga 
i us: ys iote>- 
' pace Tix 7" pone aue 3eRaY3 ; 
, s A ne 4 A’ #0. ‘MASHORG $ her a3 O27 a] ary 
OVA PRTTIRY OTE 02 ae 
| re KT Be aa 
ery WO} eH 114 
( be 
iy _ PBOaTHE 
i ~ ae ra 
* . AMOR 4A at. Fal RADSTHE, 
LA, fa ag oe Pua ere tnt ti the 
~ Ff ): £ } > s JW er 
. CCtheeeed .tes r} ova JAV TODS . 
| LUE CE VE, eH). OUR O) CBee 
Cte] | Path yu * JdG0 FED bee 
| | | Cf) At ee 
It SEP TS) POU ars) er) 
PM) COR AT Ph 2 + 3S CCR SUMED « 
ai BUD he ree i 4 JR, CCS Utd? « 
{ ¢-UG50) )ee heer 
- 2Ne8edh | 
in caine OO D- 90 
ye tey 
% <¢r Lde=(hj)U > 
a ‘te ae 
E\C T+) ee Che (1.70) (fe eOF4 . 
we Sad Lh SSuniTHOD 00% 
- wa Jo eee 2) Zoe lA 
ei) ; st 


7 v4 


| ta, tea ong 00 | 
$8.18 Aza O68 OO £ 
, BS. 49,889 pas od ie 

(00,30,.390. cos oo 
oe is, d=2 008. 
23.594387 aos 0 
er, Jn49°8. 2 & 
| oy. JH Galt 

i 00 oa 


2i,if Hel. 
7 & 40,18 
rab | 7 Jn, {,») of oq 
24.4105 208 Oo 
a= YUA py 
at woos C 


QED Cre Oe > Gace 


174 


SUBROUTINE DATARN(X,N, ERROR, NRUN) 


READS IN ONE RUN OF DATA 
CARDS, READ ACCORDING TO FORMAT STATEMENT NUMBER LAS 
CONSIST OF AN 1D NUMBER, A CARD 
FROM 1),AND TEN DATA NUMBERS, 
LAST CARD OF RUN MUST BE BLANK, 
BAST THO CARDS (OF SET MUST BE “B 
CARD OF RUN IS BLANK, OR WHEN RUNS ARE MIXED, 
WHEN: ERROR IN CARD SEQUENCE, SE 
COUNTS THE NUMBER OF POINTS IN 
REMOVE STATEMENT NUMBER 1000 IF 


Hee) 


19 


6 
1009 
qi 
8 


9 


INTEGER ERROR 
DIMENSION X(1) 
ERROR=0 
CO = 17-10 
X(J)=0.0 
CONTINUE 
BEAD GS LOT): NRUNL ET, (XC h), 1 
Paci Gor .0:) STOP 
DORIC? 70 
Peel Peel eh 20) eGe et 0. 2 
CONTINUE 
WRITE(6,1902) NRUN 
WRITECG, 106) LF,(XC1),1=1, 10 
Gila = 25.200 
L=19*J 


CONTINUE 

READ(5,191) NORUN,LTI, (X 
LECT. CTO), GOTO. 10 
N=eL 

GO) 10.9 
IFCNORUN.EQ.NRUN) GO TO 
WRITE(C6, 103) 


SEQUENCE NUMBER(CSTARTING 


LANK.. STOPS UNEEN  PRSTOAL 


TS ERROR=Wed. 
akg NN 
PR PNT OUT HOt Orsi ED 


=1,19) 


) 


CoG lisk pe) 
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WRITE(6,106) LT,(X(1),!1=K,L) 


STOP 25 
PECL oG pal f G0 210.95 


WRITE(6,104) 
ERROR=4 


WEI TEC, 106) Lh (XCl I, take) 


GO 1063 

LT=LTI 

DO 6 ako t 
LECX(1).EQ.9.9) GO 

CONTINUE 


TO 8 


WRITE(G,106) LT, (X01), 13K, 0) 


CONTINUE 
N=t-1 

WRITE(6,106) LT,(X(1),1=K,N) 
RETURN 
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ale jis) 

MOM EORMAL CT 2, LX tae OFT LO) 
102 FORMAT(1H9,'RUN NUMBER ',16//) 
103 FORMAT(1HO,10X%,'ERROR SURROUTINE DATARN = RUNS MIXED") 
104 FORMAT(1HO,10X,'ERROR SUBROUTINE DATARD - ! 
Toe "DATA CARD SEQUENCE') 
106 FORMAT(11X,18,10F10.0) 

END 
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SUBROUTINE DCLERT (CH, xX) 


PRANTS OUT N_ DATA POINTS IN _8 COLUMNS 


IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION X(1) 
N3=N/8 
L=N8*7 
WRITE(6,100) 
DO 1 t=1,N8 

K=f4+b 

silad GW SR OT a Ne Bel BB Rabe al rs We GeO A ers Fe al nes 2 
CONTINUE 
L=MON(N,8) 
NS G8 I) pal at Gai 
CONTINUE 
1=N-L+l 
DOES. v= N 

WRITE(6,102) J,X(J) 
CONTINUE 
RETURN 
FORMAT (1HO) 
EORMALGG GIA, LB! POLI ..5 3) 
EORPMALCLLOA, baad eDLL.D) 
END 
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SUBROUTINE: GRAPH: (A,B, C/DVE,N,M, IC) 


PEOTS FIVE VECTORS ON ONE GRAPH, 
EACH VECTOR HAVING TTS OWN SCALE. 
WTS: THE LENGTH OF THE. VECTORS 
MIS THEGHUMBER (OF CLINES ‘GRAPH “EXTENDS 
UNS tHe xe OPREGTLION 
LCOS (CHE (NUMBER -OF (CURVES (TO SBE IRLOTTED 
ON THECONE GRAPH 
WRITTEN BY W. VANDENBORN 


REAL *<SALB8 pOpD,e 
BRGEGER e210 PTC5), JBL, IPLCL1O) 
DVMENSHON SACL). AGB CD POO) JOCTRSECL)F ISCALECS.) 
1 BP(10) 
DATA Keens 3 ee '} ae aE ae prey Ley 
HaGlCnGhavoet se, ORePrelG. LE P.1). "GO 10 200 
1=6 
GO FO 162520015, 10,5) ,1C 
 CARPSPLOTCE, SCALE ,BOT,!1,N) 
He -GAUGARLOTCORSCALE, BOT DIN) 
TSaGrith PLOT(C, SCALE, BOT, 1,N) 
20 CREE PIUOT (8, SCALE, BOT, 1,N) 
25 GALETEPEOTCARSCALE, BOT, 1,N) 
DOTSO" te=1, LC 


IM=5-1C+] 
DOVES) J=1;,10 
h8 BP(J) = BOTCIM) + J*SCALECIM) 


Vette Ghe1000) PRTC) BOTCIM), 9 BP 
1001 FORMAT (3X,A1,1PI1E10.2,/) 
50 CONTINUE 
WRITEGGs, 1005) 
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60 
65 
70 
15 


80 


1004 


90 
1006 


100 


150 
200 
1005 


10051 


HORTOOTISSLPNOESKIR 
1=6 
DOSS bls 
PREGULY & 
GO 10 C807 
b=(-1 


BL 
Lo pig Oe OO) AC 


VE RCINT COSNGECE(S))=BOTCI))/SCALECI) *10, 


al 


OP 

[=!-1 

UPL 

l=1-1 

hPL 

IX =1X+1 

PEnakxe3G0. 5) FeGO-7F0 390 

WRITE(6,1094) [PL 

GMAT 82 SBI UNE TOA Tee Mg} 

GO TO 100 

WRITE (6,1006) IPL,J 

POR Miia (COLT axe sia” LOAT,| | Melt ei) 

fe 

CONTINUE 

WRITE(6,1903) 

RETURN 

WRITE (6,1005) 

FORMAT (' IMPROPER USE OF GRAPH- ', 
"NUMBER OF CURVES IS INVALID") 

GO TO 150 

END 


| 
VPIPCINT COSNGLCDCJU))=BOTCI))/SCALEC! ) 410. 
Gl VEGSNGL( CC) )-60T (1) )/ SCALCECI )=207, 
GIN RCUGSNGLGE CUI IRBOT CV SASECALE (1 «10. 
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SUBROUTINE PLOT, (X50 SCALE, BOT, I,N) 


PROVIDES. SCALES. FOR EACH. CURVE IN GRAPH. SUBROUTINE 
WRITTEN BY W. VANDENBORN 
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100 


REAL S. YMAK YM ERD 
DENS TORN x) wCALEK.5 7601 C5) 
ead Paull 
YMIN=X(1) 
YMAX=X(1) 
DOs aL OH 
PEC IATT (. Glas. An aL) 
Ue ede OA TD 
COC UE 
YMA =SNGLCYMAX ) 
YM! = SNGLCYMITN) 
SCALE(CI)= CYMA-YME)/19.9 
B=10.**(INTCALOGIO(SCALEC! 


) YMIN=XCJ) 
) YMAX=X(J) 


PI) 
SCALE Gl) BLOAT CIN TCOCALE CT) 7 CB? 20 ) onl) 
SEALECI OW] fell *h* OCALEN I’) 
PECSCAPEC) lr. ( YMA-YMID) (SCALEC))=SCALECI) 
1F(S5.0*SCALECI).GT.CYMA-YMI)) SCALEC! )=SCALE 
bE OEE. O00) 7G" 10° a8 
BOTCID=SCALECIO*ELOAT CUNT CYMIL/SCALECE) 9) 
60° 10 200 
BOLD = SCALE CI YELOATCINT CoM (SCALE CI) )=1) 
RETURN 
END 
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SUBROUTINE RNG(N,R,Z) 


c 
C GENERATES N PSEUDO-RANDOM VARIABLES 2 EROM A NORMAL 
C DISTRIBUTION WITH ZERO MEAN AND UNIT VARIANCE. N MUST 
C BE EVEN. R IS COMPUTED FROM TWO VARIABLES TAKEN FROM A 
C UNIFORM DISTRIBUTION (9,1) USING THE METHOD OF BOX AND 
GC PEMULLER, CANN. SOFEMATH.STAT., 12, (1958), P.610. 
C THE UNIFORM DISTRIBUTION 1S GENERATED BY A SHUFFLED 
C MULTIPLICATIVE CONGRUENTIAL METHOD. 
C MACLAREN AND MARSAGLIA,J. OF THE ASS.EOR COMP, MACH., 
Gai, (1958) 4 Be 8s 
C Z IS INITIAL INTEGER OF THE UNIFORM DISTRIBUTION 
C (DOUBLE PRECISION) 
C 
DIMENSION R(1),U(100) 
REAL G,F,Z*8 
INTEGER K,L,J,N 
C 
C GENERATE THE NUMBERS FROM A UNIFORM DISTRIBUTION 
c 
DO Ow = Wer 00 
Z=DMOD(101.0*Z+1.0,109000900.0) 
UC J) =SNGL(Z) 
10 CONTINUE 
DO W220 d=1' AN 
Z=DMOD(101. 9*Z+1.9,100000090.0) 
K=!DINT(Z) 
L=MOD(K,100090)/1000+1 
R(J)=UCL)*1. 0E-8 
U(L)=SNGL(Z) 
20 CONTINUE 
c 


C COMPUTE THE NUMBERS FROM A NORMAL DISTRIBUTION 


DO 30 Jd=1,N,2 
G=6.2831853*R(J+1) 
F=SOQRT(-2.0*ALOG(R(J))) 
R( J) =F *COS(G) 
R(J+1)=F*SIN(G) 

30 CONTINUE 

RETURN 

END 
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Typical Output of Program M2 


690000") 


910000" 0= 6 2202T*0 


VeF Cle GC 
O°O 


8T666e °O 


OCLE9C9G "OD 


C6€35506°0O 


000900000° T 


COTSE7SC°O 


SLNIUd “ON Va8v SATLVI3ea> 


8SGO9C°O = BeSOTO’O 
LeELeco* 

BOC ect 0 x -clGe Lio 
0°90 0 =C 
CILO7 76. “LCV CLEA 
B8S6T00°O = *A3d 


GOSELO"O 


9 mee BO 8) 


T9727 660-0 


Se) 


98£0GE°O 


SW BLS 


1S RN Sh Ao 7 op am Chey 2 


CO ULLve = GS 


OU GC LY Ss Ve 


OR, 


OG CSESoo "Tt 


1A do 
bLI9* SO 


CIEE AG 


0°O 


YS 


6 


A 


NOLTINSiadd L4N 


TO 


ESGESS "7G = "7 


G847l96°S 


CSetréeéey 


CC Gul oe 


0.0 


GEIVSeE.G 


"ON AVAd 


SO OV ise py 


Ch OS ESS 


SO Gaclut at 


20 


CO Ue Ca a7 


SudLinvavd 
XJONI NOLINS LA 


L wsedWwAN NN AYTAWAS 


wVad 


VLiVO 


ky 


5 


, 
’ 
ri) 
wd 


ep 
S° 


oF 
an 


. 


o ey ; 


af 


“9 343 


= 


“O° VASESL 


GO - O-geROs8e 


. “Tog aa0 


182 


S62 C2 LO" 0 


OF6901070 


E7VOLCO 40 


O50 


GETS7UC 0 


vo USBOL*Z ETEI¥BO00O 
LO GOvdEers 81S667000°0 
EO G8¥02"6 G°0 
urd 009000000°0 
70 GSEEGTS L00060000°6 
"AZG GSYWNOS 30 wWNs yguv SA1LVq34 


LLGLOOOFO 


S090£00°O 


C8col00 "0 


O°O0 


GS9¥LTO0°O 


ec6890C "Ol MLecTOeS0 ©O-G0S9S 8:9 


9S9O1TVO0 “GC 6S6110) 20 -CO-cig cle 2 


STEY COC 70 EIOLO0C. 0 cOU-dCUg ec. 1 


Oro oO 0°0 


LVL VOW Se LeSUUO "OF O=G90LCu ” | 


96£00°9 


Oe FOO *O 


02000 °O 


Cig 


NOTINALSY LAN 


cO-G 7675 Gere POO 75762)" G S 
COS-O76CSL 71 TO GvV6CEL Ss y 
7O-UGYY 9G - 9: TO OO LG 3 ac ia 

GeO Oe el 
GOs0C79S LG AO 200 2S 9 = te T 


SNOTLVIAIG DYVONVEsS 


YALIWVUVd 


| -@eeorog.! 


A aa0o% Ou 


a sTOsaconr. 2. EL1oToRO.6 


*" 
’ 


-. £0 081093.5 


ere if2.0 


183 


Typical Output of Program M3 
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TABLE VII 
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TABLE VIII 
Parameters of Power-Density Spectra 


Used to Determine the Intensity “of Jitter® 


With Thermal Conductivity Detector as Source 
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(mv) (uv?) (sec) (uv*/Hz) (uv? /Hz) 
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With Mercury Cell as Source 


aid: 1800 242 4.3 2.4 
E20 83 240 bee a6 
16.8 64000 241 24. 33% 
27.0 89000 LFS aD. 114. 
33 / 53000 240 ELS. a ee 


“craphs of the power-density spectra of both sources show 


that the spectra were roughly proportional to the Signal mean at low 
(This is indicated also by the increase of power with 


increasing mean.) Jitter and dead-time variations account for only a 
small fraction of the increase at low frequencies. Therefore, the 
increase at low frequencies, which is where the noise is most intense, 


must be due to the integrator. 


frequencies. 


peewee density at 1/2 and 1/4 in reduced frequency units. 
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